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SUMMARY 


The NASCAP code dynamically simulates the charging of 
an object in a specified plasma environment. It is fully 
three-dimensional, and it can solve complex problems in a 
few hours of computer time or less. The current contract 
called for extension, validation, and application of NASCAP. 

Numerous extensions were made in the code. They fall 
into three categories: a greater range of definable objects, 

a more sophisticated computational model, and simplified code 
structure and usage. The bulk of this report documents these 
extensions. 

An important validation of NASCAP was performed using 
a new two-dimensional computer code (TWOD) . Also, an inter- 
active code (MATCHG) was written to compare material para- 
meter inputs with charging results. 

The first major application of NASCAP was performed 
on the SCATHA satellite. A detailed shadowing study and a 
charging calculation were completed. NASCAP was installed 
at the Air Force Geophysics Laboratory, where researchers 
plan to use it to interpret SCATHA data. 
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INTRODUCTION 


This is the final report on work performed by Systems, 
Science and Software on Contract NAS3-21050, "Extension, Vali- 
dation and Application of the NA.SCAP Code" . The work was per- 
formed between September 9, 1977 and January 11, 1979. 

Most of the material contained in this final report 
was originally produced for monthly progress reports. Some 
detailed documentation has been added, and some papers pro- 
duced separately have been included. Additional documents 
produced under this contract include a revised NASCAP User's 
Manual (SSS-R-78-3739 (DRAFT) , NASA CR-159417) , and a SCATHA 
Experiment Shadowing Study (SSS-R-78-3658 (DRAFT) ) . Other 
documents of interest regarding NASCAP are: 

"Three-Dimensional Dynamic Study of Electrostatic 

Charging in Materials", Interim Report, SSS-R-78- 

3124. 

"A Three-Dimensional Dynamic Study of Electrostatic 

Charging in Materials", NASA CR-135256. 

"NASCAP User's Manual", NASA CR-13’5259. 

The above publications show the development of NASCAP and give 
background information which is not included in this report 
or in the NASCAP User's Manual. A summary of current NASCAP 
capabilities is provided in Chapter 2. 

While the first version of NASCAP, developed under 
Contract NA.S3-20119, was largely successful, it was apparent 
that many shortcomings had to be overcome to make NASCAP a 
truly useful engineering and scientific tool. These short- 
comings fell into the general categories of (1) reliability 
and ease of use; (2) generality, particularly of object 
definition; and (3) facilities for study of scientific ex- 
periments and other charging-related phenomena. 
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The use of NASCAP was made simpler and more flexible 
through expanded use of keyword input and user-specified pro- 
gram logic- These changes are described briefly in Chapter 3, 
and more fully in the NASCAP User's Manual (CR-159417). The 
NASCAP potential solver was simplified and made more reliable 
by implementation of element-by-element residual summation 
and a Scaled-Conjugate-Gradient iterative scheme. The ad- 
vantages of these techniques are discussed in Chapter 4 . 

The most difficult problem solved during this contract period 
was the failure of NASCAP' s explicit timestepping procedure. 

The LONGTIMESTEP feature (Chapter 5) was developed to guarantee 
reasonable results even when the various physical processes 
had widely disparate time constants. 

NASCAP object definition was extended in several ways. 
Thin booms (Chapter 6) extending beyond the inner mesh were 
incorporated to facilitate modeling of satellites having long 
appendages. Thin plates (Chapter 18) now provide improved 
modeling of large solar panels. Cell subdivision (Chapter 7) 
was implemented to improve resolution of an object surface. 
Also, the "patch" building blocks (Chapter 12) were added 
to simplify definition of objects having complex surface 
patterns . 

Code generality was further enhanced by taking account 
of two additional physical processes. Space charge due to 
the ambient plasma (Chapter 10) may be included in a Debye 
screening approximation. In case of materials having sub- 
stantial surface conductivity (Chapter 11) , the effects of 
this property can be evaluated. 

In addition to improving the simplicity, reliability, 
and generality of the charging simulation, facilities were 
added to study the consequences of charging. The DETECTOR 
feature (Chapter 8) allows detailed study of particles inci- 
dent upon a surface cell. The related EMITTER feature 
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(Chapter 9) may be used to study the consequences of charged 
particle emitters. A discharge analysis (Chapter 12) may be 
used to detect discharge sites and predict the effects of 
these discharges on spacecraft potential. The SHEATH option 
(Chapter 12) invokes a first-order calculation of the space 
charge density due to emitted low-energy electrons. 

Four further code development tasks were performed 
under this contract. Chapter 14 describes the conversion 
of NASCAP for a CDC 6600 computer, and its installation at 
Air Force Geophysics Laboratory. Chapter 13 describes the 
interactive MATCHG code. MATCHG treats backscatter and ' 
secondary emission from materials in a manner identical to 
NASCAP, and is intended for preliminary assessment of mate- 
rial properties. A preliminary version of a two-dimensional 
(R-6) spacecraft charging code capable of accurately predicting 
currents and space charge in the photosheath was developed for 
the purpose of comparison with NASCAP. Chapter 17 describes 
this code and the results of its NASCAP comparison. Finally, 
the HIDCEL routines were developed into a fast, highly ac- 
curate shadowing code (Chapter 13) used to produce the SCATHA 
Experiment Shadowing Study. 

Another major effort undertaken for this contract was 
development of a model of the SCATHA spacecraft for use in 
NASCAP and performance of a SCATHA charging study. The 
SCATHA application is described in Chapter 16 . 

At the close of the contract period, NASCAP was deemed 
to be in a form suitable for general distribution. A work- 
shop was held at NASA/Lewis Research Center, December 12-14, 
1978, attended by representatives of government and industry. 
This workshop was designed to introduce the attendees to 
NASCAP’ s methods and capabilities, and provide them with 
hands-on experience in its use. NASCAP is being made avail- 
able through COSMIC. 
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2 . NASCAP CAPABILITIES 


Chapter 2 of this report is a verbatim reproduction of 
a paper given at the USAF/NASA Spacecraft Charging Technology 
Conference, 31 October 1978. This paper gives a good summary 
of the form of NASCAP as it exists at the end of the contract 
year. 
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THE CAPABILITIES OF THE NASA CHARGING ANALYZER PROGRAM* 


I. Katz, J. J. Cassidy, M. J. Mandell, 

G. W. Schnuelle, P- G. Steen 
Systems, Science and Software 

J . C . Roche 

NASA- Lewis Research Center 
ABSTRACT 

Desirable features in a spacecraft modeling code are 
enumerated. The NASCAP ( NASA Charging Analyzer Program) is 
discussed in terms of its approach to the problem. Samples of 
problem setup and output are provided which demonstrate the 
ease with which the program can be used. A simple but inter- 
esting case of spacecraft charging is examined and other ap- 
plications are discussed. 

1 . INTRODUCTION 

The basic concerns of a computer spacecraft model can 
be broken down into five areas. 

1. Features of the spacecraft itself. 

2. Features of the environment. 

3. The spacecraft-environment interaction. 

4. Man-hours to set up and computer time to run a 
calculation. 

5 . A way to verify the model . 

In modeling the spacecraft itself, the point is to get 

in as much detail as can reasonably be included. This will 

vary depending on the type of model being used. The features 

desired ^ are first, some geometrical detail, such as the 

basic shape of the spacecraft body and any protrusions such 

as booms and antennae. Second, one would want to include 

which parts of the surface are bare conductor and which are 

dielectric coated. Third, it would be nice to have some 

* ‘ *1 
This work supported by the National Aeronautics and Space Ad- 
ministration, Lewis Research Center, under Contract NAS3-21050. 
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representation of the electrical circuitry connecting parts of 
the spacecraft surface. 

It is also important to decide what approximations go 
into the environment surrounding the spacecraft. The most 
basic decision is how to model the ambient plasma. Can you 
include the region far from the spacecraft, and get a detailed 
look at the region close in? Can you specify normal and ex- 
treme conditions? Does the plasma change in time? Other as- 
pects of the environment that are of concern are the sun, the 
plasma sheath, and particle trajectories. 

The spacecraft-environment interaction is mainly a matter 
of particle currents to and from the spacecraft surface. The 
important charging currents are 

•1. Incident electrons 

2 . Photocurrent 

3. Incident protons 

4. Secondary electrons from electron impact 

5. Secondary electrons from proton impact 

6. Electron backscatter 

These processes vary around the spacecraft surface, depending 
on local potential, surface material, and solar illumination. 

An ideal model would take all this local information into con- 
sideration when calculating particle fluxes. 

Computer time for spacecraft modeling can be prohibitive. 
A model that is general ends up solving a series^ of equations 
with hundreds or thousand’s of variables . An exact solution is 
enormously expensive, and it may be hard to get convergence 
from an iterative solution. Much care must be 'put into this 
aspect of the problem, lest an otherwise elegant modeling pro- 
gram start to impersonate an infinite loop. 

The most expensive way to verify a modeling program is 
to build a spacecraft like the model and send it up. Other, 
more reasonable techniques, are to model ground experiments. 
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to check answers for reasonableness, and to test the program 
on known problems. 


2 . NASCAP APPROACH 

As we have seen, the physics which must be examined in 
order to model spacecraft charging presents a problem of for- 
midable dimensions. It would be impractical to develop a com- 
puter code that was state of the art in every aspect of the 
problem. By placing restrictions on the class of problems to 
be examined we have been able to construct the NASA Charging 
Analyzer Program which provides useful information in those 
cases of most practical interest. It is most applicable to the 
high voltage charging caused by magnetospheric substorms. 

Our approach has been to limit the range of ambient en- 
vironments to those whose Debye lengths, A^, are large compared 
to object dimensions. For magnetospheric substorms this is 
definitely true. 

8 e m 10,000 eV 

n . -3 

n ml cm 
e 

A d m 0.7 km 

Only for the very largest conceivable spacecraft are object 
dimensions comparable to Debye lengths. For finite Debye 
lengths we have included ambient plasma screening approxima- 
tions, albeit of modest applicability. 

Overall, we have modeled all aspects of the problem ex- 
cept electromagnetic wave propagation. Our idea has been to 
use the best available analytical theories wherever possible 
and to minimize the brute force number crunching . By doing 
this we have been able to combine good treatments of ambient 
environment, sheath, complex object, and electrical and parti- 
cle interactions into a single code. This is done by using 
known physics and developing approximate models where necessary. 
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For example, NASCAP contains analytical approximations to 
electron backscatter as a function of electron energy and angle. 
While not as accurate as Monte Carlo transport results , these 
formulations do give reasonable yield estimates and can be 
evaluated quickly at hundreds of surface locations each time- 
step. Thus we obtain reasonable estimates in reasonable amounts 
of time as opposed to best estimates regardless of cost. This 
philosophy permeates the code. Where quasi-analytical models 
were necessary but unavailable, we have developed them. 

The procedure followed in the code is to approximate the 
spacecraft in a 3-D Cartesian grid. Free space around the satel- 
lite is provided by nesting grids within grids where each grid 
has a linear dimension twice that of the grid it surrounds. 

There can be an arbitrary number of these nested grids. How- 
ever, the more grids, the longer the computer time per calcula- 
tion (see Figure 1) . 

All parts of the spacecraft must remain in the innermost 
grid, except for booms which can extend into several grids. 

The object itself is composed of an assembly of cubes, sliced 
cubes, plane surfaces, and skinny cylinders, as shown in Fig- 
ure 2. Each surface can be of an independently specified mate- 
rial, with up to 15 different materials permitted (Figure 3) . 
Certain classes of surfaces may be subdivided for higher reso- 
lution. 

Object definition is by far the most complicated aspect 
of using a three-dimensional computer code. To make the pro- 
gram easy to use, NASCAP provides an extremely simple object 
definition language. Complex three-dimensional spacecraft can 
be described with a minimum of effort. The satellite shown in 
Figure 4 is a good example. The central structure is octagonal 
with a gold circumference and aluminum top and bottom surfaces. 
The two planar sheets represent solar cells with kapton cover- 
ing the back surface. They are attached to the main body with 
kapton coated cylinders. This object was defined using 31 brief 
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lines of input (Figure 5) . The simple object definition com- 

[21 

mands are fully explained in the NASCAP User's Manual. 

Once the object definition is complete, the program al- 
ternately calculates charge accumulations on surfaces, and 
potentials caused by these charges. Due to the variety of time- 
scales in the system, the algorithm used to advance the charge 
distribution in time is extremely complex, so complex that it 
uses a couple thousand element self-generated capacitor model 
as its own internal estimator. 

NASCAP produces a variety of printed and graphical out- 
put. The fundamental idea is to help the user follow the 
progress of the calculation (Figures 6-14) . 

The first graphic output is a two-dimensional view of 
the spacecraft with surface cells shaded to show the material 
types. Each surface cell is individually classified by mate- 
rial, with up to 15 different material types allowed. 

Next is a three-dimensional perspective view of the 
spacecraft without hidden line removal. This is helpful in 
tracking down object definition problems. It is followed by 
a view from the same perspective with surface cells outlined. 

In this surface cell plot, hidden lines are removed. The user 
gets a quick and accurate feeling for the defined object. The 
routine that generates these plots also calculates exposed sur- 
face areas for determining photoelectron emission. 

These plots are generated at object definition time, be- 
fore the actual satellite charging begins. The major outputs 
of the charging calculation are the flux breakdown printout 
and potential contours. 

The flux breakdown printout shows, for any surface cell(s) , 
the charging currents operating on that cell. Each individual 
surface cell requires a separate calculation. By requesting 
flux breakdown printouts, the user can closely follow the 
charging process at any point on the surface. 
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Contour plots are an efficient way to show what’s hap- 
pening to the electrostatic potential both near the spacecraft 
and far away. The user can look at the potential contour plots 
generated every time cycle and get a good feeling for global 
changes in the spacecraft sheath. 

NASCAP detector routines plot flux density versus energy 
of particles reaching the detectors. Detectors can be placed, 
at the user's discretion, on any surface cell. 

The emitter routines plot trajectories of particles 
emitted at various energies. These trajectories, along with 
potential contour plots, give a very good idea of fields sur- 
rounding the spacecraft or test tank object. 

Finally, if local electric field stresses exceed some 
user specified threshold value, a message is printed and the code 
redistributes charge as if a discharge had occurred. 

3 . VALIDITY OF THE MODEL 

With a model as broad in scope and as complex (over 400 
subroutines) as NASCAP, the immediate question is "How do you 
know that it gets reasonable answers?" So that we have confi- 
dence in NASCAP results, testing and comparing to analytical 
results has been a major part of the development program. The 
accuracy of the various components have been examined in con- 
figurations simple enough to determine their inherent accuracy. 

Since the capacitance of simple objects such as spheres, 
cubes and cylinders are known quite well, we have used these to 
determine how well the potential routines work. For all cases 
the NASCAP results were within 10 percent of analytical pre- 
dictions, and for objects of more than a zone resolution and 
for booms of radius much less than the grid spacing, the NASCAP 
results were accurate to a few percent. The electric fields in 
space were of corresponding accuracy near the satellite and in- 
creasing accuracy away from the vehicle. The accuracy of the 
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potentials are limited only by the ability of the finite element 
interpolation functions to represent the true solution. For 
complex objects, the NASCAP code uses the same algorithms and 
the accuracy should be comparable. Since NASCAP automatically 
takes into account mutual capacitances, it is a vast improve- 
ment over hand generated capacitor models for complex spacecraft. 

NASCAP assumes that charge is accumulated on, as opposed 
to deposited within, dielectrics. Bulk conduction is included. 

We have performed detailed one-dimensional calculations of charge 
transport within dielectrics, and have found this to be a reason- 
able approximation for electrons of a few to tens of kilovolts in 
all but the thinnest of dielectrics. It is also an approximation 
that can easily be modified in the future if the need arises. 

The charging currents are the algebraic sum of incident 
fluxes and backscattered, secondary, and photoemitted electrons. 
For spherical test cases we have compared NASCAP reverse tra- 

[31 

jectory currents with spherical probe formulas. Depending 

on the number of trajectories sampled the results were in 
reasonable agreement, the largest errors due to the differences 
between numerical and analytical integrals over angle of the 
backscatter and secondary emission formulas. Thus the two basic 
requirements of a charging calculation, the potential and charge 
accumulation, are performed well by NASCAP. 

The NASCAP material interaction models have been developed 
from literature results. Their predictions are being compared 
with laboratory experiments and are the subject of another paper. 
It should be pointed out, however, that NASCAP accepts para- 
meters for these models as input and that the models themselves 
are contained in very short, easily replaceable subroutines. 
Consequently, modifications and improvements in the formulations 
can be made very simply if needed. 

The particle trajectory algorithms are second order ac- 
curate in particle timesteps insuring good conservation of 
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energy and magnetic moment. Orbits are followed beyond the 
outermost grid boundaries by using an extrapolation of the mono- 
pole potential. This allows long excursions of emitted parti- 
cles to see if they return to the spacecraft. 

The algorithm employed to integrate charging currents 
over a timestep is quite complex to ensure physical results. 
Rather than describe the technique in detail, we present a cal- 
culation which illustrates how it works. 

A simple example, which nevertheless displays some of 
NASCAP ' s usefulness as a model, is the case of a spherical ob- 
ject in sunlight. Since the photocurrent is larger than the 
incident electron current, a capacitor-current balance model 
would lead one to the conclusion that a sunlit surface will re- 
main at a positive potential relative to the surrounding plasma. 
However, the NASCAP charging current integration routines recog- 
nize that space charge limiting prevents photoelectrons and 
secondary electrons from supporting a potential barrier of more 
than a few volts. This feature, combined with the multidimen- 
sional aspects of the potential leads to a very different equi- 
librium, one with the illuminated surfaces a kilovolt negative. 

We ran NASCAP for the case of a teflon coated sphere in 
sunlight. The environment for this case is an isotropic, 
Maxwellian plasma with a temperature of 20 keV and a density 
n e = n^ = 1 cm ^ . Sunlight was incident on one side of the 
sphere (Figure 15) . 

Figures 16-22 show the time development of the electro- 
static field. (The satellite-sun line lies in the plane of 
these figures. Dark and sunlit cells are differentiated by 
shading.) For the first a.0.1 second the sphere charged uni- 
formly. Over the next few seconds, the negative charge accumu- 
lated by the shaded surfaces began to dominate the electrostatic 
field, causing a saddle point to appear in front of a sunlit 
surface. At about 10 seconds the potential at the saddle point 
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became negative. The sunlit surface maintained a potential a 
few volts positive relative to the saddle point. Final steady 
state is reached with the sunlit surface at -1.0 kV and the 
shaded surface at -3.6 kV'. 

The final steady state potentials were reached at time 
4 

t ^ 10 sec. This involved some 30 timesteps, and used total 
computer time of about one-half hour. Thus in a reasonable 
amount of computer time NASCAP can provide good physical in- 
sight into charging phenomena, insight which is unobtainable 
using simpler computer models. 

4. APPLICATIONS OF NASCAP 

NASCAP is designed primarily to give engineering esti- 
mates of spacecraft potentials during magnetospheric substorms. 
It also can provide detailed particle spectra for a given en- 
vironment and spacecraft potential configuration in order to 
aid in interpreting results of scientific experiments. As of 
this time the applications of NASCAP have been limited to the 
comparison with laboratory material charging test results and 
to the generation of models of a few scientific spacecraft. 
Comparisons have been done to validate the material properties 
portion of the code. 

One application of NASCAP which is of engineering im- 
portance is the study of active charging control. The opera- 
tion of onboard charged particle beams has been proposed as a 
means of minimizing the effects of ambient environment space- 
craft charging. NASCAP features an emitter algorithm that 
models the trajectories and charge transfer effects of such 
beams. For example, we have placed a one kilovolt, one milli- 
ampere electron emitter on a satellite precharged to -2.5 kV. 

The potentials on spacecraft ground and on an insulated surface 
as a function of time are shown on Figure 23. Notice that the 
insulator will differentially charge to a substantial negative 
potential. Sample particle trajectory plots during the charging 
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phase are shown in Figure 24. By modeling such systems NASCAP 
can estimate their utility and point out any severe design 
problems, so that actual flight experiments have the best 
chance for success. 

An important problem, particularly in the future, is the 
interactions of large space structures. While not specifically 
designed for this application, the finite Debye length sheath 
treatment in the NASCAP code will combine with the reverse tra- 
jectory particle flux routines to give good estimates of space 
charge limited charge collection. The present algorithm employs 
linear Debye shielding (Figures 25-26) . In the future, models 
of the ambient plasma sheath more relevant to dense collision- 
less plasmas, will be implemented. The object definition rou- 
tines can already handle objects of large size by decreasing 
the object resolution (Figure 27) . 

The most ambitious application to date is the generation 
of the SCATHA model. This model utilizes the full capabilities 
of the code. The model and some preliminary calculations are 
the subject of another paper. 
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A two-dimensional view of the first four nested meshes. Each succeeding 
mesh increases the volume of calculation space by a factor of eight. Ca 
culation time is roughly linear with the number of meshes. 





Figure 2. NASCAP can simulate virtually any object that can be 
built from these fundamental shapes — cube, three 
types of sliced cube, planar square, and thin cy- 
linder. 
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.Figure 3. The spacecraft surface is made up of as many as 1200 
surface cells. Each cell is assigned a material type 
and an underlying conductor. The surface cell may 
represent either bare conductor or dielectric layer. 
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Figure 5. Object definition. The object in the preceding 
figure (paddle satellite) is defined by these 
commands . 
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Figure 7. Object structural plots give a- perspective view 
without hidden line removal. 
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Figure 8 



Surface cell hidden line plots give a clear idea of 
overall spacecraft structure. 
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Figure 9. A breakdown of charging currents can be requested for any surface cell. 
This information is given at each timestep. 



Figure IQ. Two-dimensional potential contour plots give a 

clear picture of electrostatic potential at each 
timestep. 
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Figure 11. Particle emitters can be specified at any surface 

cell. This plot shows particles from five emitters 
for various angles of emission. 
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Figure 12. Particle detector plots show energy versus flux 
density. Detectors can also be located at any 
surface cell. 
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Figure 13. Graphic output for the test tank case includes trajectories of electrons 
from the source to the object. 
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Figure 14. Potential contours around a fully charged teflon covered grounded plate in 
a ground test tank. An electron beam is coming from the left. Notice the 
fully formed potential saddle point to the right of the plate. 




Figure 15. A NASCAP sphere — modeled as a tvzenty-six faceted 
object. This one is 3 meters in diameter with 
158 surface cells and 144 surface nodes. 
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Figure 16 . 


Potentials on shadowed and solar illuminated sur 


faces of a teflon sphere in a plasma (N = lO^/m 
6 = 20 keV) . 









Figure 18- Potential contours around sunlit sphere showing 
early appearance of saddle point (x) at -5.6 
volts . 
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Figure 19. Potential contours around sunlit sphere showing 

fully formed saddle point at approximately -8 volts. 
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Figure 20 . Potential contours about sunlit sphere showing 
saddle point at approximately -25 volts. 
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Figure 21. Steady state potential contours about sunlit sphere. 
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POTENTIAL (KILOVOLTS) 



Figure 23. Active control simulation. A 1 rnA particle emitter 
is activated with beam energy of 1 keV. The space' 
craft goes from a negative 2.5 kV potential to posi' 
tive 1.0 kV. Spacecraft ground remains at about 
that level while a solar cell on the surface falls 
back to a negative potential. f 



Figure 24. Particle emitter trajectory plot. Some of the 

emitted particles escape the spacecraft vicinity, 
while others return to various points on the sur- 
face . 
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Figure 25. An approximate screening expression is employed to 
show shielding effects. Shown is a two meter cube 
charged to -100 V, in a plasma with Debye length 
of 33 meters. 
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Figure 26. Here the same cube is charged once again to -100 V. 

This plasma has Debye length of 3.3 meters . The 
denser plasma leads to more significant shielding, 
and the potential falloff is steeper near the cube. 
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3 . INPUT 


In March of 1978 the NASCAP input routines were sub- 
stantially rewritten, resulting in greatly simplified running 
procedures. The changes were, first, modular program control, 
and second, default keyword values. 

The total input required from a typical NASCAP user 
is first, a set of command words, and second, a set of three 
input files which describe the object, the environment, and 
the user-selected program options. 

The command words allow modular program control. The 
user can perform object definition or not, or ask for a shadow- 
ing calculation or not, before any charging analysis is per- 
formed. User options can be changed between charging timesteps. 
The control word names and their functions are given in the 
NASCAP User's Manual (DRAFT) . ^ 

The files used to describe the object and the environ- 
ment remain essentially unchanged. The user options file is 
new. All of the program control options which used to be in- 
cluded in the NASCAP runstream are now in the options file. 

These are quantities like file numbers, graphics and printed 
output control, length of timesteps, and size of computational 
grid. Default values are supplied for all of these options. 
Options not specified assume the default values. This im- 
provement has greatly simplified NASCAP use. Options and 
defaults are listed in the NASCAP User's Manual (DRAFT) . ^ 
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4 . POTENTIAL SOLVER 

Two major changes were made in the NASCAP potential 
solver. The first was to calculate the coproduct Au in terms 
of the volume elements rather than in terms of the nodal 
points. As a result, coding for 11 special elements" is greatly 
simplified. Boom-type elements and other new types can more 
easily be included. The second change improved the Poisson 
solver routines. The old conjugate gradient technique has 
been replaced by a scaled conjugate gradient technique. The 
new method takes approximately the same amount of computer 
time per iteration, but converges in far fewer iterations. 

4.1 ELEMENT BY ELEMENT COPRODUCT 

The element oriented coproduct (residual) calculation 
required major restructuring of many routines in the TRILIN 
section of the code. It was accomplished in the following 
manner. 

The major task was to form the residuals for each 
potential value from volume or surface "stiffness matrices" 
operating on the potentials, as opposed to combining the ele- 
ment stiffness matrices together into a giant matrix. The 
advantage is that the local matrix bandwidth is reduced 
greatly. Additionally, NASCAP 1 s I/O time has been greatly 
reduced, with only a modest increase in CPU time. 

Algebraically, NASCAP now generates the residual 
vector r by the product 

r = Ad) , 
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where 


A = £ Dj , 

elements {j} 

Uj being the stiffness matrix for element j. connects only 

the nodal points that bound element j , and for a simple element 
is a symmetric 8x8 matrix. A is a sparse matrix with typi- 
cally 27 nonzero entries in each row. Since only the residual 
vector r is used in the conjugate gradient potential solver, 
we need not store A but rather form the residual vector from 
element residuals, r ^ . 



Algebraically, the two techniques are identical 

J = X) 5j = X (u ji> = ( X u j) 4 = • 

j j x j ' 

But operationally within NASCAP, the formation of residuals 
element by element greatly simplifies the coding necessary to 
treat booms, struts, thin plates and surface subdivision. 

4.2 SCALED CONJUGATE GRADIENT 

The changes necessary to implement the scaled conjugate 
gradient method were mainly in the POTENT (subset of TRILIN) 
section of the code. This method improves convergence tre- 
mendously in cases of large zone size and very thin dielectric 
skins. -It is accomplished by scaling the large coproduct 
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matrix so all diagonal elements are of order unity. We de- 
scribe below three conjugate gradient methods for the solution 
of linear equations Mx = y: (1) the original NASCAP method; 

(2) a scaled method; and (3) a simple computational method equiv- 
alent to (2) . Method (3) is used in the present version of NASCAP. 


Method 1: The Ordinary Conjugate Gradient Method. 
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Method 2: The Scaled Problem. 
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Finally, 


n rWT>"l \ n 

x - D (D x) 



Notice that in this code formulation the solution vector x and 
the matrix M need never be explicitly scaled. 


48 



5. LONGTIMESTEP OPTION - SUBROUTINE LIMCEL 


During the past several months the analysis performed by 
the LIMCEL routine, invoked by specifying the LONGTIMESTEP option, 
has become a cornerstone of the NASCAP modeling effort. Through 
this analysis the applicability of NASCAP has been extended to 
physical regimes where we had not previously dared to attempt 
calculations; such as differential charging in sunlight and 
active control by high current emitters. In this chapter we 
discuss the need for the LIMCEL routine and the logic used in it, 
as well as some of the physical principles upon which the logic 
is based and some of the mathematics used to apply it. 


5.1 NEED FOR THE LIMCEL ROUTINE 

At the end of our first contract year a serious short- 
coming of NASCAP was apparent: NASCAP could not satisfactorily 

simulate charging of an object initially dominated by emission 
of low-energy secondary- or photo-electrons. Other problems 
were also seen from time to time, such as a tendency for poten- 
tials to oscillate unstably in time and space, and an inability 
to handle materials with substantial conductivity. 

The root of these problems lies in the disparate scales 
of time and distance (i.e., capacitance) which typify the charg- 
ing process. The capacitance per unit area of a satellite is 
given by 



s 

z 10 pf/m 


2 


(5.1) 


where A is the satellite area and R its effective radius. Thus, 

-5 2 

with a typical charging current of 10 A/m a satellite will 
charge at a rate 


V = J/C v 10 volts/sec. 


(5.2) 
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However, the process of true interest is differential charging, 
characterized by the thickness of dielectric coatings, d ^ 10 - ^ m, 
leading to 




0.1 lif/m^ 


(5.3) 


Thus a similar current level produces 
V D = J/C D 'v 10 0 volts/sec. 


. (5.4) 


From Eg. (5.4) we see the desirability of performing simulations 
on a 0.1 - 10.0 second timescale, while Eq. (5.2) says that on 
such a scale small changes in net current will lead to wide, 
non-physical oscillations. 


5.2 AN IMPLICIT CHARGING TREATMENT 

The considerations of the previous section forced us to 
abandon the original "explicit" charging treatment used in the 
first version of NASCAP in favor of a more stable "implicit" 
algorithm. In simplified form, the basic equations are: 

Explicit: C[V(t 2 ) - V(t^)] = J(t^)(t 2 ~t^) (5.5a) 

Implicit: C[V(t 2 ) - V(t 1 )] =J(t 2 )(t 2 “t) (5.5b) 


The obstacle to solving Eq. (5.5b) is that, while J(t^) is known, 
J(t 2 ) is a complicated function of the unknown, V(t 2 ). If, how- 
ever, we make the approximation that 


J(t 2 ) = J(t x ) 


+ J' (V(t 2 ) - V{t 1 ) ) 


Eq. (,5.5b) gives 


v(t 2 ) - v(t x ) 


J(t x ) (t 2 -t x ) 

C - J' (t 2 -t x ) 


(5.6) 


(5.7) 
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If we take a case of C = 10 ^ f, J(t n ) = 10 ^ A, t_ - t.. = 

_ g J- Z i. 

1 sec, J' = -10 A/volt (J' < 0 is required for physical sta- 
bility) we find 

Explicit: V(t 2 ) - V (t^) = 10 5 volts (5.8a) 

Implicit: V(t 2 ) - V(t^) = ^9.9 volts (5.8b) 

That (5.8a) is unstable while (5.8b) is stable is indicated by 

-3 

plugging into (5.6), giving J(t 2 ) =10 A (explicit), or 
J(t 2 ) = 10 ^ A (implicit) . 

Implementation of the implicit algorithm in NASCAP is 
made more complex than solution of scalar equation (5.5b) by 
the matrix-vector nature of the charging problem. Far more dif- 
ficulties, however, are raised by careful assessment of Eq. (5.6). 
For some processes, such as linear surface and bulk conductivity, 
(5.6) is exact and the matrix J' is known. For others, such as 
current collection from a plasma, it is adequate at best, and 
determination of J' is not a trivial matter. For still others, 
such as cutoff of low energy emitted electrons or emitter cur- 
rents by potential barriers, (5.6) is totally inadequate. When 
these processes dominate, we must take the approach of estimat- 
ing the final potential based on known conditions and deter- 
mining a mean current consistent with the final potential. 

5.3 OVERALL STRUCTURE OF THE LIMCEL ANALYSIS 

The objective of the LIMCEL analysis is to determine ap- 
propriate coefficients and boundary conditions for, and to solve, 
the equation 

C [ V(t 2 } - V(t x )] = [J(t 2 ) + 0V(t 2 )] (t 2 -t 1 ) (5.9a) 

which can also be written (using (5.6)) 
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[V(t 2 ) - Y (t l )] = ? (t l } + ? v(t l ) (5.9b) 

Here J denotes the net current excluding current due to conducti- 
vity processes. On entry to LIMCEL the capacitance and conducti- 
vity matrices C and o are known, as is the initial potential 
V(t^) ■ The explicit external current J(t^) is also known, al- 
though some recentering is required. Additional information 
used by LIMCEL includes potentials at points external to the 
satellite, the portion of J(t 1 ) due to low energy emitted 
electrons, and information concerning active control emitters 
on the spacecraft. After an optional discharge analysis, LIMCEL 
returns the left hand side of Eq. (5.9a) to TRILIN for use in 
updating the total charge vector and POTENT then performs the 
full calculation of new potentials on and about the spacecraft. 
Within LIMCEL, Eq. (5.9b) is solved (usually several times) 
using the "Incomplete Cholesky Conjugate Gradient" (ICCG) 
method. ^ 

The LIMCEL analysis then proceeds in the following six 
phases (see block diagram. Figure 5.1): 

1. Preliminary Phase . Prior to beginning the timestepping 
process a lumped-circuit-element model (Figure 5.2) of the 
spacecraft is constructed. The nodes of this model are the 
conducting satellite segments (maximum of 7) and the grid points 
and subdivision points located on insulating surfaces (maximum 
of 1024) . (A further restriction is that the number of non-zero 
matrix elements in Eq. (5.9b) may not exceed 9537.) 

2. Explicit Phase . The right hand side of (5.9b) is 
evaluated using known information. Also, data concerning low 
energy electron emission, external electric fields, and emit- 
ters is processed for later use. If no LONGTIMESTEP analysis 
was requested, control is returned to TRILIN. 
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Figure 5.1. Block diagram of LONGTIMESTEP option. Subroutine LIMCEL is enclosed by the 
dashed line. 
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3. Preliminary Charging Analysis . Those nodes for which ap- 
proximation (5.6) is likely to be in error are identified and 
preliminary estimates of their final potentials are made. A 
set of "trial" potentials for all the nodes is formed for use 

in the next phase. 

4. Flux Derivative Determination . The fluxes to the surface 
cells are calculated using "trial" potentials for the purpose of 
determining the matrix J'. The flux derivative matrix is assumed 
diagonal . 

5. Final Charging Analysis . Equation (5.9b) is solved 
repeatedly for the final potentials, and the left hand side 

of (5.9a) is evaluated for the mean currents. Those constraints 
found to be unnecessary are removed. Nodes which remain con- 
strained are set at potentials consistent with their mean cur- 
rents. The left hand side of (5.9a) is evaluated for return 
to TRILIN. 

6. Discharge Analysis (Optional) . The discharge analysis 
is described elsewhere. 

Throughout the LIMCEL segment use is made of the poten- 
tial limiting input parameter, DVLIM. The "trial" potentials 
and flux derivatives are found consistent with the notion that 
DVLIM is the maximum potential change desired. If, during the 
final charging analysis, any conductor displays a greater po- 
tential change, the LIMCEL analysis is repeated with a shortened 
time step. 

5.4 PRELIMINARY PHASE: LUMPED CIRCUIT MODEL 

The spacecraft lumped-circuit model constructed by 
NASCAP is shown schematically in Figure 5.2. The nodes of the 
circuit represent either conducting segments of the spacecraft, 
or points located on dielectric surfaces. The circuit elements 
are : 
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Figure 5.2. Lumped-circuit model of spacecraft constructed by 
NASCAP and used for LONGTIMESTEP option. 


w > f 


a. "Small" capacitances, coupling circuit nodes to 
plasma ground. 

b. "Large" capacitances, coupling surface nodes to 
conductors, or conductors to each other'. 

c. Resistors representing bulk and surface conducti- 
vity processes. 

The matrix [C - cr] corresponding to this network is sparse, 
symmetric, and positive definite — ideal for treatment by 
the ICCG inversion algorithm. 

Formation of the circuit matrix begins with subroutine 
GENMTL, called automatically following OBJDEF from NASCAP 
through DRISCM. GENMTL first forms the PTLIST array (see 
Figure 5.3) listing all grid points and subdivide points on 
dielectric surfaces. (Boom nodes are added later by' B00M2 . ) 
The PTLIST establishes the index numbers of the circuit nodes 
in the circuit potential vector; conductor nodes are indexed 
sequentially following the surface nodes. Next the "matrix, 
skeleton" is formed. The skeleton is an integer array indi- 
cating the non-zero matrix elements: the negative -entry (-i) 

corresponds to the ith diagonal element, and is followed by 
an arbitrary number of positive entries (j) (in order) indi- 
cating off-diagonal matrix elements between nodes i and j . 
GENMTL forms only those rows of the matrix corresponding to 
PTLIST nodes; rows corresponding to conductor nodes are added 
later by LSTMAT, taking into account conductors which are 
biased or held at fixed potential. Finally, GENMTL forms 
the surface conductivity matrix, a s * 

Further development of the circuit model occurs in 
response to the CAPACI keyword. Subroutine CAPACI calls 
POTENT to calculate the potential about the spacecraft with 
a unit charge on spacecraft ground. This information is then 
used by CELGET to calculate the overall body capacitance, C a , 
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Figure 5.3. Surface node list (PTLIST) entry format. 
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and to apportion it among the circuit nodes: 


C 

CO 



8 * 

= — / 

v J 


■+ A 

E • n dS 





(5.10) 


where the sum runs over surface cells. The surface cell quan- 
tities are recentered among the' surface nodes and conductors 
to obtain the circuit matrix quantities. Before exiting CELGET, 
the far right of Eq. (5.10) is renormalized to the far left. 

The renormalization factor is printed, and .should be unity to 
within about 1 percent. (A renormalization factor substantially 
different from unity should not be accepted if the user plans to 
use the LONGTIMESTEP option.) Subroutine CELGET also forms the 
array of "large" capacitances for the surface nodes, and the 
bulk conductivity matrix elements. 

Included in the PTLIST array are surface nodes having 
capacitive/resistive coupling to more than one conductor. These 
multiconductor points (maximum of 128) are indicated by a zero 
conductor index in their PTLIST entries. These entries are 
duplicated in the array MULTCN, which serves as a map to the 
arrays CMULT and SIGMLT in which the matrix elements are stored. 


5.5 EXPLICIT PHASE 

The explicit phase of LIMCEL (subroutines ADEMIT, GETDQ, 
and DQSCND) is concerned with transforming the relevant informa- 
tion available in NASCAP to arrays paralleling the circuit node 
list. ADEMIT is concerned with modification of charging by 
active control particle emitters. In cases where particles re- 
turn to the satellite surface the fluxes are appropriately 
modified. The height and location of potential barriers seen 
by emitted particles, together with other information, is 
stored in the /EMITR2/ common block. 
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GETDQ is the major routine of this phase. It forms the 
following arrays: 


VPTS 

DQO 

DQEMIT 

EPTS 


the potential on each surface node. 

the net explicit charge accumulation on each cir- 
cuit node, including emitter current and bulk con- 
ductivity. 

the low-energy electron current emitted from each 
circuit node. 

the effective electric field at each node: 

EPTS { j ) = s(i) E. • n ± / £s(i) (5.11) 

i i 

where the sum is over surface cells, and s(i) is 
the low energy electron current from cell i at- 
tributable to node j . 


Finally, DQSCND modifies the explicit charges DQO to reflect 
the explicit contribution of surface conductivity. 


5.6 PRELIMINARY CHARGING ANALYSIS 

The preliminary charging analysis serves two functions: 
(1) to identify those circuit nodes for which Eq. (5.6) is an 
invalid approximation ("bad” nodes), and (2) to determine a set 
of "trial" potentials for use in the next phase. The first 
operation is to find the "bad" nodes. A node is marked "bad" 
if 

1. It is a conductor to which a particle emitter is 
grounded, or 

2 . Its explicit current is caused to be positive by 
low energy electron emission, and its potential o£ 
effective electric field is positive (electron 
attracting) . 

Those nodes which have not been marked bad but are nonetheless 
electron attracting have their low energy emitted electron cur- 
rent reduced. Trial potentials are determined for the "bad" 
nodes such that secondary electrons or high energy emitted 
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currents will marginally escape, subject to the constraint 
that the potential change can be no greater than DVLIM. The 
ICCG potential solver is then called, and the preliminary 
trial potentials inspected for additional "bad" nodes. (For 
example, an uncharged, sunlit satellite initially has no "bad" 
nodes. However, it may return from ICCG charged many kilo- 
volts positive, so that all photoemitting nodes will then be 
marked "bad".) If additional nodes are found "bad", ICCG 
generates a new set of preliminary trial potentials. Finally, 
the trial potentials are formed by subjecting the ICCG result 
to the constraint that no node may have a potential change 
larger than DVLIM. 


5.7 DETERMINATION OF TRIAL POTENTIALS FOR FIXED NODES 

When a circuit node is marked "bad" its potential is 
to be fixed such that the offending particle can marginally 
escape. Specializing to electrons (for ions, sign change 
and interchange of min and max are appropriate) , that means 


V. 

i 


V . + 

mxn 


e/e 


(5.12) 


where e is an energy characterizing the particle emission, and - 
V min t ^ ie potential occurring on an escaping parti- 

cle trajectory at the next times tep. During the preliminary 
charging analysis the energy e is set to zero for low-energy 
electrons and near the maximum emitter energy for emitters . 

This provisional value is to be refined during the final 
charging analysis. 

Determination of the potential V . presents a more 

mxn 

•difficult problem. Since V . occurs somewhere in space, it 
c mxn 

is clearly beyond the scope of the circuit model. We ap- 
proach this problem by dividing it into two questions: 


1 . 


What is 
V ( ? ) ? 


the current value of 


mxn 


this potential, 
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2. How is this value likely to change in proceeding 
to the next timestep? 

In answering these questions the emitter case is clearly the 

easier, since particle tracking information is available to 

give the value of and the location of sheds light 

on the second question. For the low energy emission case we 
guess 


V ( ? } = min 
mm 


0, 


v! 0) 


- E< 0) 


Ax 


(5.13) 


where is the current node potential, the current 
effective field [Eq. (5.11)], and Ax the mesh spacing. We 
then address the second question by supposing 


V . 
min 


v (°) 

mxn 



so that Eq. (5.12) becomes 


V. 

i 


V ^ 
mm 






+ e/e 


(5.14) 


(5.15a) 


or 


V. 

1 


v <°> 

mm 



a) . 


(5.15b) 


For the low energy case we choose a = 0 for V m ^ n = 0 , and 

a = 1/3 for V . <0. (The value 1/3 is usually an underesti- 
' mm 

mate, but a £ 1/2 tends to produce instabilities.) For the 
emitter case we choose 


a = 0 . 7 e min ( 1 , R_,/R 0 ) 

U -D 

where R is the spacecraft's capacitive radius, R- the radius 

C / n \ C 

at which V. was found, and 0.7 has been inserted to insure 
mm 

stability. Finally, the potential is subject to the constraint 
|V. - V. (0) I < DVLIM . 

1 l l 1 


61 



The approximations involved in the foregoing treatment 
tend to produce a "multidimensional lag" in the timestepping 
process. For example, the appearance of a saddle point at 
timestep 8 will not be reflected in the surface potential until 
timestep 9 , at which point it may be substantially underesti- 
mated. One gains, however, the advantage of being able to 
proceed stably toward a steady state with minimal use of the 
expensive POTENT routine and a minimal amount of particle 
tracking. One or two short timesteps are usually sufficient 
to resolve the multidimensional lag. 


5.8 FLUX DERIVATIVE DETERMINATION 

The purpose of including J' in Eq. (5.9b) is to assure 
that those properties of nature which cause physical stability 
in real satellites also provide mathematical stability in 
NASCAP. As was the case in the previous section, uncertainties 
are resolved in favor of providing additional mathematical 
stability. 

The coding, then, proceeds as follows: 

1. The trial circuit node potentials are used to find 
trial surface cell potentials (subroutine VSHARE) . 

2. Incident, backscattered, and secondary fluxes are 
calculated using the trial potentials (subroutine 
GETFLX) . 

3. The trial fluxes are compared with the original fluxes 
passed from TRILIN (subroutine GDQTRI) . 

4. The difference between the trial and original fluxes 
are recentered to the circuit nodes (subroutine 
QSHARE) . 

5. The current differences are divided by potential dif- 
ferences to obtain provisional values for the diagonal 
matrix J’ (subroutine DFDV) . 
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6. Final values of J* are determined to assure stability 

(subroutine DFDVMX) : 

For surface nodes: 

J* = min (J’, -| J(t 1 ) |/(2*DVLIM) ) 

For conductor nodes : 

J' = min (J 1 , 0) . 

5.9 FINAL CHARGING ANALYSIS 

The purpose of the final charging analysis is to deter- 
mine the charge accumulation on each circuit node consistent 
with Eq. (5.9b) and the information known concerning the "bad" 
nodes. A flow chart of this portion of LIMCEL is shown in 
Figure 5.4. 

First ICCG1 is called to solve Eq. (5.9b), and IMPFI 
evaluates the left hand side of (5.9a). At this point con- 
ductor potentials are checked, to see if a timestep reduction 
is in order. Next HIGHQ examines charge accumulation on "bad" 
surface nodes to see if any are unphysically too far positive. 

If such nodes are found, they are unfixed and the analysis 
begun anew. If not, the "bad" surface nodes are examined to 
see if any are accumulating more negative charge than would 
be the case for total cutoff of low energy electron emission 
(subroutine LOWO) . If such nodes are found, J(t-^) is re- 
placed by the value it would have in event of total cutoff 
of low energy emission, the node is unfixed, and we once again 
begin the final charging analysis. 

Thus when subroutine VFIX is executed all "bad" surface 
nodes satisfy 

J(t x ) > J > J(t 1 ) - L (5.16) 
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where L is the low energy electron emission current and J the 
mean current calculated by IMPFI. VFIX determines the energy 
e of Eq. (5-12) such that 

J = J(t x ) - L ( 1 - e~ z/2 ) (5.17) 

and redefines the constrained potentials of the "bad” nodes 
accordingly. The "2" appearing in (5.17) corresponds to the 
2 eV characteristic energy of photo- and secondary-electrons. 
Subroutine VCFIX performs a similar function for conductor 
nodes, except that a more complex formulation than (5.17) is 
used for conductors to which emitters are grounded. 

Since the constraint potentials have now been readjusted, 
ICCG1 and IMPFI must be called once again,- and various unfixing 
checks made (UNFIX, UNFIXC) . If more nodes are unfixed we 
once more again restart the final charging analysis, but this 
branch is taken only once. Otherwise, we are done. 

Output of this section are the predicted potentials of 
all the circuit nodes and the corresponding charge accumulation. 
The predicted potentials are used in the optional discharge 
analysis. If discharges are found, the charges are modified 
accordingly. The charge accumulations are then returned to 
TRILIN for use by subroutine QUPDAT . 

5.10 THE DECREASE IN EFFECTIVE PHOTOCURRENTS DUE TO SADDLE 

POINTS IN ELECTROSTATIC POTENTIALS NEAR DIFFERENTIALLY 

CHARGED SPACECRAFT 

This section is a verbatim reproduction of a paper given 
at the 1978 IEEE Annual Conference on Nuclear and Space Radia- 
tion Effects, Albuquerque, New Mexico, July 18-21, 1978. The 
authors are M. J. Mandell, I. Katz, G. W. Schnuelle and P. G. 
Steen, Systems, Science and Software, and J. C. Roche, NASA- 
Lewis Research Center. 
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THE DECREASE IN' EFFECTIVE PHOTOCURRENTS DUE TO 
SADDLE POINTS IN ELECTROSTATIC POTENTIALS 
NEAR DIFFERENTIALLY CHARGED SPACECRAFT* 

I . INTRODUCTION 

As interest in spacecraft charging has grown over the 
past decade, many spacecraft charging calculations have ap- 
peared in the literature. Such calculations may be character- 
ized, roughly in order of increasing complexity, as 

1. Equilibrium current balance calculations.^ ^ 

Ml 

2. One-dimensional computer programs. 

3. Lumped-circuit-element computer programs. 

r g — 8 1 

4. Multidimensional computer programs. 

The first type of calculation simply predicts the floating 
potentials of surfaces having particular material properties in 
various environments. Such calculations demonstrate that space- 
craft can indeed charge to high negative potentials, and deter- 
mine the relative importance of various material and environ- 

M1 

mental properties. One-dimensional codes introduce the ad- 
ditional complication of a photoelectron sheath which can sub- 
stantially modify the dynamics of charging and the final poten- 
tial distribution. 

r 5 1 

Lumped-circuit-element codes model a complex satellite 
electrically as a network of capacitors and resistors. By 
assigning to each node a current-voltage characteristic M (V^) , 
a dynamic charging calculation can be performed. However, since 
the code has no geometrical knowledge of the satellite, effects 
such as shadowing, incomplete particle trajectories, particle 
reflection, and photosheath effects are either totally neglected 
or inserted "by hand” . 

ye 

This work was supported by the National Aeronautics and Space 
Administration, Lewis Research Center, Cleveland, Ohio under 
Contract NAS3-21050. 
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The purpose of this paper is to illustrate the presence 

of important multidimensional effects in spacecraft charging. 

r s 1 

Two-dimensional codes have been under development by Parker 1 

f 7 1 

and by Laframboise . The calculation described below was 

performed using the three-dimensional NASA Charging Analyzer 
Program (NASCAP) . 1 J NASCA? dynamically simulates the charging 
of an object made of conducting segments which may be entirely 
or partially covered with thin dielectric films. The object 
may be subject to either ground test (electron gun) or space 
(magnetospheric) environments. The simulation alternately 

(1) treats the accumulation and emission of charge by surface 
materials and its redistribution by conduction processes, and 

(2) calculates the electrostatic potentials on the object and 

in the surrounding space. Implicit algorithms allow simulations 
of long periods of time, and particle tracking capabilities 
enable calculation of such quantities as response of charged 
particle detectors. NASCAP also has extensive graphics capa- 
bilities . 
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II. STATEMENT OF PROBLEM 


A. TEST OBJECT 

The "satellite" to be charged is a 3-meter diameter 

-4 

sphere whose surface consists of a 10 m teflon coating over 
a conducting substrate. In NASCAP, a "sphere" is modeled as 
an object having 26 faces. This calculation was performed on 
a "sphere" having 158 surface cells and 144 surface nodes 
(Figure 1) . 

Some of the material properties ascribed no teflon are 
given in Table I . It is worth noting that the conductivity 
value, which is larger than indicated by low field measurements, 
may be appropriate to the equilibrium electric field of "ilO 
volts/meter . 

B. ENVIRONMENT 

The environment: was an isotropic, Maxwellian plasma 

appropriate to a severe magnetospheric substorm, having a tem- 

-3 

perature of 20 keV and a density n^ = n^ = 1 cm . This plasma 
has a Debye length of one kilometer, so that the space charge 
contribution of the ambient particles is totally negligible. 
NASCAP was run in a mode in which each surface cell collected 
incident fluxes of electrons and protons appropriate to a 
spherical probe at the local surface potential. Sunlight was 
incident on one side of the sphere. 
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Figure 1. NASCAP representation of sphere used in this cal- 
culation, showing surface resolution. 
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TABLE I 

TEFLON PROPERTIES USED IN THIS CALCULATION 


Dielectric constant 
Thickness 

Conductivity (bulk) 
Conductivity (surface) 
Effective atomic number 
Effective atomic weight 
Density 

Secondary yield-electron 
impact 

6 

max 

E 

max 

Secondary yield-proton 
impact 

Yield for 1 keV proton 

Energy for -maximum yield 

Photoemission (normally inci- 
dent sunlight) 


2.0 

—d 

10 ' meters 

-1 d 

10 ‘ (ohm-m) 

(neglected) 

10 

16.7 

2 . 2 gm-cm 2 
3.0 

0 . 3 keV 
1.4 

70 keV 

2 x 10 _o A/m 2 



III. APPROXIMATE PHOTOSHEATH MODEL FOR STRONG 
DIFFERENTIAL CHARGING 


The sheath of low-energy electrons which can form near a 
positively charged surface is known to have complex structure, 
dynamics, and transport properties. NASCAP has the capability 
of determining photosheath currents through tracking of emitted 
particles. However, not only is such a procedure time-consuming, 
bun it jeopardizes the numerical stability of the calculation. 
This is because photocurrents are sensitive to surface potential 
changes comparable to the two-volt characteristic energy of 
emitted electrons , and thus small compared with the kilovolt 
differential potentials of interest. The purpose of this section 
is to justify a principle which can be used to determine the 
potential of photoemitting surfaces. To this end, we first 
show that any substantial electric field can dominate space 
charge effects in determining photosheath structure. It then 
follows that the surface potential will attain a value such that 
the fraction of photoelectrons escaping over an electrostatic 
barrier is just that needed to maintain current balance. 

Let us then consider space charge-limited emission in the 

presence of an external field. If the field is negative (i.e., 

into the surface) no sheath will form, so we will treat only 

positive fields. For the simple case of monoenergetic (energy E) 

electrons emitted normally from a plane surface, a virtual 

cathode will form at a distance d from the surface. The sheath 

[ 9 ] 

thickness d is found using the space charge equation 

/ dv \ 2 = ( dv ) 2 _ 8J / mV 
\dxj \dx) e Q \2|e 
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with the boundary condition 


V (d) = 0 
V (o) = E/je| 

= external field. 

x=d 

Fiaure 2 shows the sheath thickness as a function of external 

" 2 . 
field for the parameters J = 3 nA/cm and E — 2 eV . It: is ap 

parent that any substantial positive external field will com- 
pletely dominate space charge effects and suppress emission of 
low energy electrons . Taking into account the distributed 
spectrum of low— energy emitted electrons , we are led to the 
following principle: 

UndcA. condition* o ^ *tA.ong dl$ i CA.cntlal akaA.gd.ng a. 
pko to emitting *aA.{ s acc will A.cach a pozcntlal Mick a* to main- 
tain a positive cx.tcA.nal clcctA.lc filcld o& a &ew volt* pcA. 
mctcA- . 
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100.0 


1000.0 


Electric Field (volts /meter) 


Figure 2 


Electron sheath thickness outside a planar surface 
emitting 2 eV electrons as a function of electric 
field in the low current limit (dashed line) and 
for 3 nA/cm^ emitted current (solid curve) . 



IV. RESULTS 


NASCAP was run to calculate the electrostatic potentials 
on the surface of, and in the space surrounding, a sunlit 
teflon-coated sphere. Currents to the sunlit surfaces were 
determined based on the principle put forth in the previous 
section. From an initial uncharged state, the sphere reached 
a final steady state having 2.5 kV of differential charging. 
Figure 3 shows the potentials on a shaded and a sunlit surface 
cell as a function of time. 

Figures 4-8 show the time development of the electro- 
static field. (The satellite-sun line lies in the plane of 
these figures. Dark and sunlit cells are differentiated by 
shading.) For the first ^0.1 seconds the sphere charged uni- 
formly (Figure 4) . Over the next few seconds, the negative 
charge accumulated by the shaded surfaces began to dominate the 
electrostatic field, causing a saddle point.to appear in front 
of a sunlit surface (Figures 5-6) . Ac about 10 seconds the 
potential at che saddle point became negative. In accordance 
with the principle put forth in the previous section, the sunlit 
surface maintained a potential a few volts positive relacive to 
the saddle point (Figure 7) . Final steady state (Figure 8) is 
reached wich the sunlit surface at -1.0 kV and che shaded sur- 
face at -3.6 kV. 

The components of incident and emitted current are shown 
in Table II. It is apparent thac low energy emitted electron 
currents are always dominant on the sunlit surface. In the 
final steady scate, the net negative current incident upon the 
shaded side is balanced by conduction through the teflon, which 

7 

has an internal field of 1.0 x 10 volts/meter. 
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Potential (volts) 



Figure 3. Potentials on shadowed and solar illuminated sur 
faces of a teflon sphere in a plasma (N = 10^/m 
6 = 20 keV) - 







Figure 4. Potential contours about a sunlit sphere early in 
time . 
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Figure 6 . Potential contours around sunlit sphere showing 

fully formed saddle point an approximately -8 voles. 
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Figure 7 


Potential contours about sunlit sphere showing saddle 
poinr at approximately -25 volts. 
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Figure 8. 


Steady state potential contours about sunlit sphere. 
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Source 


TABLE II 

2 

COMPONENTS OF INCIDENT AND EMITTED CURRENT {A/m ) 

Shaded Side 


Sunlit Side 


Incident Electrons 

Resulting Backscatter 


Incident Protons 


Photocurrent (L) 


Total 


Initial 

Final 

Initial 

Final 

-3.78 x lO" 6 

-3.16 x 10“ 6 

-3.78 x 10 -6 

-3.58 x 10 -6 

1.02 x 10" 6 

8.49 x 10~ 7 

1.02 x 10~ 6 

9.62 x lO" 7 

1.23 x 10 -6 

1.03 x 10" 6 

1.23 x 10" 6 

1.17 x 10 -6 

8.83 x 10 -8 

1.04 x 10” 7 

8.83 x 10~ 8 

9.31 x 10~ 8 

9.15 x lO" 7 

1.09 x 10 -6 

9.15 x 10 -7 

9.66 x 10~ 7 

0 

0 

2.00 x 10 -5 

2.00 x 10'" 5 

-5.34 x 10“ 7 

-9.45 x 10 -8 

1.95 x 10~ 5 

1.96, x 10” 5 


(L) - Indicates low energy emitted electrons subject to space charge or field limiting 




V. DISCUSSION 


It is instructive to consider a lumped-circuit element 
solution to this problem (Figure 9) . The shaded and sunlit 
surfaces are coupled to spacecraft ground by resistance R and 
capacitance C^, and to plasma ground by capacitance C p . Since 

^ s R^/d and C„ v £ R where a is che dielectric thickness 
Do P o 4 

and R the satellite dimension, C Q ^ 10' C p . The current to the 
sunlit surface is dominated by che emitted photoelectrons, 
which are absent from the shaded side- As indicated in the 
figure, for R = the surfaces at equilibrium maintain their 
individual floating potentials, while finite resistivity 
ameliorates somewhat the degree of differencial charging, but 
leads to no qualitative differences. It is only when multi- 
dimensional effects, manifested through the electric fields 
external to the satellite, are taken into account that a sunlit 
surface can develop a negative potential. 

To further illustrate the saddle point effect, the 
problem was rerun under the assumption that the low-energy 
electrons are emitted with a characteristic energy of 50 volts. 
(Actual photoelectrons have energies of about 2 eV. However, 
the 50 volt choice leads to results more suitable to NASCAP 
spatial resolution.) The final potentials were -700 volts on 
the sunlit side and -3400 volts on the shaded side. Potential 
contours for this case are shown in Figure 10. The saddle point 
can be seen more clearly than in Figure 8 because of the higher 
positive electric field outside the sunlit surfaces. Tra- 
jectories for electrons of one to one hundred eV energies are 
shown in Figure 11. It is apparent that all the low-energy 
electrons return to the surface, while the highest energy 
particles escape over t-he saddle point. 
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Figure 9. Lumped-circuit-element model and solution (schematic) for charging of 
sunlit sphere. On the I-V plot, the crosses indicate current balance 
for R = «>, and the triangles for finite R. 


* 

Figure 10 . Steady state potential contours about a sunlit spnere 
whose low-energy emitted electrons are assumed to haye 
a characteristic energy of 50 eV. 
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Figure 11. Trajectories of electrons emitted from two different 
surface cells in potential of Figure 10. Electron 
initial energies are logarithmically spaced from 
1 to 100 eV. 





VI. CONCLUSIONS 


This relatively simple example illustrates the importance 
of multidimensional effects in spacecraft charging. In particu- 
lar, the mechanism demonstrated here may have played a role in 
events where satellites ATS-5 and ATS- 6 were observed to charge 
negatively in sunlight, whereas simple current balance would 
have predicted that photoemission would keep the sunlit side, 
if not the entire satellite, at a positive potential. Thus sim- 
ple current balance calculations can lead to erroneous con- 
clusions about equilibrium as well as dynamic charging. However, 
these errors need not be due to complex plasma-sheath-dynamic 
phenomena, but may be caused by relatively simple electrostatic 
effects . 

A further conclusion is that a dielectric patch, when 
differentially charged to high negative potential, can, through 
its effect on particle trajectories, have an influence out of 
proportion to its .area. Such effects have been suspected in 
cases of spurious particle detector response. Another possible 
effect might be to prevent escape of actively emitted electrons. 
This would annul the intended discharge of a satellite by low 
energy emitters . 

In summary, multidimensional electrostatic effects play 
an important role in spacecraft charging. Surfaces at high 
negative potential can suppress emission of low energy electrons 
elsewhere. Such effects can explain observation of negative 
potentials on sunlit surfaces, and may seriously affect particle 
detector response and active potential control mechanisms . 
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6 . BOOMS 


Booms are skinny cylinders — less than one grid unit 
in radius. They must be an integral number of grid units in 
length. Unlike other portions of NASCAP objects, booms are 
not confined to the innermost mesh. They are also unique in 
that the boom radius is a real number between zero and one. 
This allows very accurate modeling of boom volume and surface 
area. 

The most severe restriction of the boom object is that 
booms must lie along grid lines. They can only point in the 
directions of the coordinate axes, and a boom crossing a mesh 
boundary must line up with a grid line in the outer (coarser) 
mesh. 

In three-dimensional plots, booms are shown as skinny 
rectangular parallelepipeds. But within the NASCAP treatment, 
boom circumferences are truly round. All boom calculations 
are made appropriate for a curved, not a flat, surface. 

6.1 BOOM DEFINITION 

Booms are defined by giving starting and ending points, 
a radius, and a surface material. The definition routine ex- 
pects to read 5 cards per boom, in the following order: 

CARD 1 . CCODE 

FORMAT (A6) 

CCODE must contain the literal 'BOOM' 
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CARD 2 . CCODE , IXB , IYB ,12b, IGB , IXE , IYE , I ZE , IGE 
FORMAT (A6 ,4X, 815) 

CCODE must contain the literal 'AXIS' 

IXB, IYB, and IZE give the starting coordinates of the 
axis, and IGB gives the starting grid. Similarly, the 
ending information is given in IXE, IYE, IZE, and IGE. 
Any offset applies only to coordinates in grid 1. Co- 
ordinates in outer grids are assumed to have OFFSET = 
( 0 , 0 , 0 ). 

CARD 3. CCODE, RADIUS 

FORMAT (A6,4X,F10.0) 

CCODE must contain the literal 'RADIUS' 

RADIUS gives the boom radius in inner mesh units. 

CARD 4. CCODE, MATERIAL 
FORMAT (A6,4X,A6) 

CCODE must contain the literal ' SURFAC ' 

MATERIAL must be a previously defined surface code. 

CARD 5 . CCODE 

FORMAT (A6) 

CCODE must contain the literal 'ENDOBJ' 

6.2 RESTRICTIONS ON BOOMS 

There are several restrictions which apply to booms: 

1. Booms must be defined parallel to a coordinate axis. 

2. Booms must not intersect one another, and there 
must be at least one node separating parallel booms. 

3. Booms must intercept other objects only ' at (1,0,0), 

(0,1,0), (0,0,1), (-1,0,0), (0,-l,0), or (0,0,-l) 

surfaces . 

4. Booms must not pass through objects. 


c 


- V- 
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5 . 


Booms must pass from an inner grid towards outer , 
grids . 

6 . Booms may intercept grid interfaces at right angles 
to the interface; however , at least one node must 
separate each boom from grid interfaces which are 
not intercepted. 

7. Booms may not approach within one unit or less the 
bottom of a thin plate. 

6.3 BOOM CELL FLUX SUMMARY 

Flux breakdown printouts can be obtained for boom sur- 
face cells. Boom surface cells come at the end of the surface 
cell list, continuing the cell numbering scheme. 

To request a flux breakdown printout for a boom surface 
cell, first find its position on the boom cell list. Add this 
position number to the highest number for a non-boom cell — 
the last cell printed in the standard surface cell list. For 
example, to get a printout on the tenth (10) boom cell when 
there are five hundred and two (502) non-boom surface cells, 
insert in the RDOPT input file a card reading 

SURFACE CELL 512 
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7. SUBDIVISION 


NASCAP can now subdivide user-specified surface cells 
for finer potential resolution. Specified surface cells can 
be subdivided into 9, 16, or 25 nodes, depending on whether 

•k 

the user requests 1, 2, or 3 subdivisions. 

The major difficulties in this process are those of 
data management. Implementation required forty-odd new sub- 
routines and changes to fifteen old ones. 

The new code has been tested on some simple objects. 

It gives good results (see "RESULTS") . 

The sections which follow deal with various aspects of 
the implementation . 

7 . 1 RESULTS 

A test case was run on a surface with two conductors. 

The surface was four grid units (meters) square and one grid 
unit thick. One half was held at 100 V and the other at 
200 V. 

The following graphs (Figures 7.1 - 7.3) show the poten- 
tial falloff as the conductor boundary was crossed on a line 
1 cm above the surface . 

The entire minus Y surface was subdivided to the maxi- 
mum extent (NSUB = 3) so that the 4x4 grid became an effec- 
tive 17 x 17 grid. 


The three subdivision limit is a storage limitation. The 
method is general and could be used for any number. 
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Volts 






7.2 CELLS AND ELTS 

NASCAP deals with two basic geometric types — surface 
cells and volume elements. Surface cells are basically two- 
dimensional structures about one grid unit to a side. They 
define the spacecraft volume by enclosing it. The spacecraft 
is a three-dimensional object with surface cells making up 
its skin. 

Volume elements fill space. Most of them are cubes 
with edges of one grid unit. Some cubic volume elements co- 
incide with volume occupied by the spacecraft. These volume 
elements inside the spacecraft are said to be "filled". 

Filled volume elements do not figure into NASCAP calculations. 

Sometimes the spacecraft boundary cuts across a volume 
element. Such an element is divided into two pieces. One 
piece is filled — it does not exist. The remaining piece, 
the "empty" piece, becomes a non-cubic volume element — a 
tetrahedron perhaps, or a wedge, or a truncated cube. 

Every surface cell is either the border between two 
volume elements, or the border between two parts of a single 
volume element — one filled and one empty. 

In the following pages, surface cells will be referred 
to simply as "cells" . Volume elements will be called "volume 
elts" or "elts". Cells have two dimensions. Elts have three. 
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7.3 


FACES AND EDGES 


Subdivision is initiated when a user specifies a surface 
cell or set of cells which is (are) of special interest. The 
user wants a more accurate potential value for points on the 
cell than is obtained with bilinear interpolation from the 
corners. Subdivision creates a number of new nodes in between 
the regular grid nodes. These new nodes will have a role in 
all of NAS CAP ' s potential calculations . 

By far the most complicated part of the potential cal- 
culation is the coproduct A x U = AUN. Coproduct terms are 
calculated for each surface cell and each volume elt. Cells 
and elts that have subdivided nodes must be treated in a 
special way. 

The center cell in Figure 7.4 has been subdivided. It 
has two extra nodes on each edge and a group of four additional 
nodes toward the center of the face. This cell is a face- 
subdivided cell, or FACE-SD cell. Each of the four cells 
marked "E" is subdivided along one edge, the edge they share 
with the face-sd cell. Cells with one subdivided edge are 
called EDGE-SD cells. There are special coproduct routines 
for face-sd .cells and edge-sd cells. 
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Figure 7.4. Extra nodes on subdivided cell. 
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As noted in the section CELLS AND ELTS, a surface cell 
is often the border between two volume elts — one filled and 
one empty. Any subdivided cell must form such a border. The 
filled elt is ignored. The empty elt in front of the cell must 
be included in coproduct calculations. 

An empty elt that lies against a face-sd cell is a 
FACE-SD elt. An elt lying against an edge-sd cell is an EDGE-SD 
elt. There are special coproduct routines for face-sd elts and 
edge-sd elts. 

Subdivided cells and elts have by definition extra nodes. 
Extra nodes that lie on a cell edge are called EDGE nodes. Ex- 
tra nodes not along an edge, i.e., nodes not shared by two cells, 
are called FACE nodes. 

If N is the number of subdivisions, a face-sd cell has 

2 

N edge nodes along each edge. The cell has N face nodes and 

2 , 2 

four grid nodes. The total is N + 4N + 4 or (N + 2) nodes 
on this surface. 

In summary, the extra nodes used by NASCAP for subdivi- 
sion are classified as face nodes or edge nodes. Surface cells 
and volume elements are changed by the introduction of these 
new nodes. The four new types of cells and elts are called: 

Face-sd cells. 

Edge-sd cells. 

Face-sd elts . 

Edge-sd elts . 
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7.4 


SDINPU — OBJECT DEFINITION 


The user specifies which cells are to be face-sd. 

This choice is limited by several restrictions (see "RESTRIC- 
TIONS" section) . The object definition word "FINER" initiates 
this process by causing a call to subroutine SDINPU. The 
first three integers on the "FINER" card are the origin of a 
rectangular group of surface cells which are to be face-sd. 

The next three integers are the AX, AY, and AZ of the rec- 
tangle, one of which must be zero. The final integer is NSUB, 
the number of subdivisions in each direction. NSUB is either 
1, 2, or 3, and is the same for all subdivided surfaces. If 
different NSUB's are specified, the one on the last FINER 
card is the one that counts . 

The input is in the usual 15 format starting in col- 
umn 1 1 . 

Example : 

FINER 8882012 

subdivides the surface cell with corners (8,8,8), (9,8,8), 

(8.8.9) , (9,8,9) and the one next to it (9,8,8), (10,8,8), 

(9.8.9) , (10,8,9). 
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7.5 


RESTRICTIONS 


The basic restrictions on the use of subdivision are 
listed below. A set of examples which clarify these restric- 
tions follows . 

1. Only square surface cells and empty cubic volume 
elts may be subdivided. 

2. An edge-sd cell or elt may have only one edge 
subdivided. 

3. A face-sd cell or elt may have only one face and 

. the four edges of that face subdivided. 

4. Only cells and elts in the innermost mesh may be 
subdivided, and a subdivided volume elt may not 
touch the edge of the innermost mesh. 



SURFACE CELL INCONSISTENCIES 



NOT ALLOWED 


The above configuration of face-sd cells leaves two cells 
with multiple sd edges. The next drawing clarifies the 
situation. 
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NOT ALLOWED 

Same Situation (NSUB = 2) 


The cells marked "E" are acceptable edge-sd cells. The cells 
marked "X" are edge-sd cells with two subdivided edges. A 
simple solution to this problem is shown in the next drawing. 
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ACCEPTABLE 


The cells marked "S" have been specified by the user 
as face-sd cells. The program will classify the cells marked 
"E" as edge-sd cells. In many situations, if nearby cells 
cannot be selected as single face-sd cells, the solution is 
to subdivide a larger area including the cells of interest. 
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NOT ALLOWED 


ACCEPTABLE 


The cell between the 
two face-sd cells, is 
double edge-sd. 









































MOT ALLOWED NOT ALLOWED 

Edge-sd cell tri- 
angular not square 



ACCEPTABLE 














NOT ALLOWED 

' Cells marked "X" are 
double edge-sd 





NOT ALLOWED 

Four faces of unit 
cube are double 
edge-sd 
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Be careful at corners. The reader should verify 
that if the object being considered is a unit cube in space, 
there must be 0, 1, or 6 subdivided surface cells. There is 
no legal way to subdivide 2, 3, 4, or 5 of the surfaces. 
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VOLUME EL? INCONSISTENCIES 

It is possible for the surface cells to be subdivided 
in a consistent way, while at the same time one or more volume 
elts are inconsistent. Some examples: 



NOT ALLOWED NOT ALLOWED 

A double face-sd volume elt The same elt is double 

edge- sd 



S, S'' NOT ALLOWED 
S, S" NOT ALLOWED 
S',S" NOT ALLOWED 



S, S" NOT ALLOWED 
S^,S"" NOT ALLOWED 
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There is no easy solution to inconsistency problems 
on concave surfaces. 
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SURFACE CELLS THAT ARE NOT SQUARE 

It is easy to remember that surface cells that are 
not square cannot be subdivided. But remember also that 
volume elts which have any non-square face cannot be sub- 
divided. 



NOT ALLOWED 

Edge-sd elt -has diagonal 
on one face 
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7.6 


• SDLIST 


User input selects which surface cells are to be 
face-sd. It remains to identify face-sd elts, edge-sd cells, 
and edge-sd elts. 

Face-sd elts are easily found. There are two volume 
elts sandwiching each face-sd surface cell. One of them must 
be a filled volume elt. The other must be empty. The empty 
elt is designated face-sd. 

To identify edge-sd cells and elts, a list of sub- 
divided edges must be created. There are four sd edges for 
each sd face, but some edges are shared by two sd faces. 
Redundancies are eliminated. 

Once the sd edges are listed each surface cell is 
tested to see if it has an edge on this list. If it does, 
and if it is not a face-sd cell, it is designated an edge-sd 
cell. 

Any edge in space is the intersection of four volume 
elts. For each subdivided edge, a test is run on the four 
volume elts it touches. 

If the elt is empty and not already subdivided, it be- 
comes an edge-sd elt. If it has already been noted as sub- 
divided, a check is run for consistency. An edge-sd elt mus 
have only one subdivided edge, and a face-sd elt must have 
one subdivided face and four subdivided edges . 

In addition to identifying subdivided cells and elts, 
SDLIST creates four lists which are used by many other sub- 
division routines. They are: 

An index of sd faces (LOFACE) . 

An index of sd edges (LOEDGE) . 

A list of sd cells (JSUBBR) . 

A list of sd elts (LTABBR) . 
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7.7 


RESIDUALS AND SUBDIVISION 


The most basic operation in the coproduct section is 
residual calculation. 

We know how to calculate residuals for the vertices 
of a rectangular parallelepiped with a node at each vertex. 
We do not- know how to calculate residuals for a rectangular 
parallelepiped with extra nodes lying along an edge or on a 
face. 

Therefore, we break sd cells and elts into smaller 
pieces that are easy to handle. The manner of breakup is 
shown in the following diagrams. 


Ill 



A xace^sd elt showing all node 
points (NSUB = 2) . 


Values are interpolated for 
imaginary nodes. 


So the elt can be broken into 
9 rectangular parallelepipeds 
for residual calculation. 




An edge-sd elt with all nodes 
shown QISU3 = 3) . 




Imaginary edge nodes are 
interpolated . 



For breakup into four rectangular 
parallelepipeds . 
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DICER is the subroutine 
■and elts. SLICE takes care of 


that handles face-sd cells 
edge-sd cells and elts. 
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7.8 


STORAGE 


The student of NASCAP will realize that at the time 
when COPROD calculates AUN for volume elements of the inner 
grid, there is some effort expended to make room for three 
large (17 x 17 x 33) arrays. U, AUN, and the element table 
are all necessary for this calculation. 

The AUN calculation for a subdivided element requires 
these three arrays plus SDAUN plus SDYOU plus the edge and 
face index arrays. These additional' arrays total only about 
five thousand words, but there is no place to put them. 

The solution is to create an abbreviated element table 
including only subdivided elts. Then after the usual calcula- 
tion has been made for non-subdivided elts, the element table 
is dispensed with. The space it occupied is written over with 
the abbreviated element table and the other subdivision infor- 
mation. The large U and AUN arrays remain where they were. 

The IOBJ file, formerly not used, now holds informa- 
tion for subdivision. The arrays, in the order stored, showing 
maximum dimensions, are: 

SDARR (18*MXFACE) 

SDYOU (18*MXFACE) 

SDDIV (18*MXFACE) 

SDAUN ( 18*MXFACE ) 

SDPEE ( 18*MXFACE) 

SDROUS ( 18 *MXFACE ) 

LOFACE (MXFACE) 

LOEDGE ( 3* MXFACE) 

LTABBR (2, 6* MXFACE) 

JSU3BR (4*MXFACE) 

In the current code MXFACE = 50. This is the maximum 
number of subdivided faces allowed. The implicit limit on 
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number of subdivided edges is =3*MXFACE, maximum number of 
subdivided elts =6*MXFACE, and maximum number of subdivided 
cells =4*MXFACE . 

The first six sd arrays correspond to existing CG 
arrays. LOFACE and LOEDGE are the index arrays for faces and 
edges. LTABBR and JSUBBR are abbreviated element table and 
surface cell list. 

JSUBBR is just selected entries from the JSURF list. 

The surface cell location is encoded. However, the spacial 
location of a volume element listed in the element table is 
known by its position in the element table. So for our con- 
densed version we need two words for each elt. 

Word 1 = SQUISH — encoding of elt origin. 

Word 2 = element table entry for subdivided elt. 

7 . 9 INTERNAL REPRESENTATION 

Faces and Edges 

Sd faces and sd edges are represented by single word 
entries in the sd face index list and sd edge index list. 

Any face on the .grid is determined uniquely by a normal 
direction and an origin. The normal direction is given as 1, 
2, or 3, corresponding to plus or minus X, plus or minus Y, or 
plus or minus Z direction, respectively. The origin is the 
lowest X, Y, and Z coordinate on the face. 

Any edge on the grid is similarly determined by its 
direction and origin. Notice that while this scheme uniquely 
identifies a face among faces or an edge among edges, con- 
fusion will result if faces and edges are mixed. The face 
with direction (normal direction) 3 at origin (8,8,8) is a 
square parallel to the X-Y plane with corners at (8,8,8), 
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(9,8,8), (8,9,8), and (9,9,8). The edge with direction 3 at 

(’8,8,8) is a line segment from (8,8,8) to (8,8,9). 

For subdivided faces and edges, the index word is 

formed by subroutine SQUISH. SQUISH forms an integer by 

. 6.4 2 

adding 10 * direction +10 * X coordinate +10 • Y coordi- 

nate + 10 u • Z coordinate. The reverse transformation is ac- 
complished by EXPAND. 

Cells 

The surface cell list is encoded as in Figure 7.5. The 
lowest order 5 bits hold the material number for a surface. 

Material numbers 11 - 20 inclusive' are used for face-sd 
cells. Numbers 21 - 30 are for edge-sd cells'. 

The material number mod 10 is used to reference mate- 
rial properties stored in the MATPR array. 

Elts 

The element table for volume elements is described in 
Figure 7.6. A face-sd elt is set to element type 10. An 
edge-sd element is element type 11. 
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5 43 210 987654 321098 765432 109876 5 43210 


U I l 



HGF E D C B A 


Field Bits 

A 4-0 Material index 

B 5 Set for right-triangular 100 sur- 

faces and for 111 surfaces whose 
enclosing volume cell is mostly 
empty 

C 11-6 Direction of surface normal (in 

crystallographic notation) 

D 17-12 Z-coordinate 

E 23-18 Y-coordinate 

F 29-24 X-coordinate 

G 32-30 Conductor index 

H 34-33 Orientation code for right- 

triangular 100 surfaces 


Figure 7.5. Surface cell list (JSURF) entry format. 
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CODE FOR ELEMENT TABLE [LTBL (NX,NY,NZ) ] . 


54321 0987654321 098 765 432109876 5 43210 

I III! Ill 


E DC B A 


Field 

A 

B 

C 


D 

E 


Bits 

4-0 Elt-type code 

14-6 Orientation code 

18 Set if elt is completely filled 

(interior) 

19 Set for an empty special elt 

30-21 Index used to reference PHOJ array 
to determine low energy electron 
currents 


ORIENTATION CODE 

3 X-3 bits. Each group of 3 contains 1, 2 
lower 2 bits, with the high bit set for negative. 


m m 9 m 

Code (-) i i 1 , (-) ±2 , (“) i 3 


takes (r 1 ,r 2 ,r 3 ) to {-) , (-) r^ , (-) 


e.g., the following codes take a point to (x,y,z): 


or 3 in the 



3 


Octal Code 

123 

365 

532 

176 

567 

617 


Figure 7.6. Element 


Dec ■ Code 

1,2,3 
3, -2,-1 
-1,3,2 
1,-3, -2 
-1,-2, -3 
-2, 1,-3 


table codes and 


Point 

(x,y, z) 

(-z ,-y ,x) 

( -x , z , y ) 
(x,z ,-y) 
{-x , -y ,-z) 
(y,-x,-z) 


orientation codes . 
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7.10 SIGNIFICANT SUBDIVISION CONSTANTS 


Most of these constants are computed in SDLIST. The 
exception is NSUB which is user-specified. 

NSUB 

The number of subdivision points that lie between 
adjacent grid points on a subdivided surface. 


NUMFA, NUMED 

NUMFA is the total number of subdivided faces. NUMED 
is the total number of subdivided edges (NUMED <_ 

4 • NUMFA) . 

NSDELT , NSDURF 

NSDELT is the number of sd volume elts, NSDURF the 
number of sd surface cells. 


LENTRU 

This is the number of extra nodes, i.e., the number of 

2 

subdivided points. It is the sum NSUB • NUMED + NSUB 
• • NUMFA. 
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7.11 SIGNIFICANT SUBDIVISION ARRAYS 


The first four arrays are constructed by SDLIST. They 
keep track of what goes where. The last six arrays are where 
what goes. They are the data arrays used by TRILIN, and they 
are the object of the complex indexing schemes that subdivision 
routines deal with. Unless noted, the arrays are one-dimen- 
sional. 

LOFACE, LOBDGE 

These are the index arrays used for storage of sd 
point information. Each sd face or sd edge has. a 
one-word entry in the table. The entry is an encod- 
ing {from SQUISH) of the low index grid node of the 
face or edge, and the direction of the face or edge. 

The direction of an edge is the direction it points — 
the direction of a face its normal direction. 

From the location of an entry corresponding to the 
sd face or edge, the storage location of its points 
(in SDYOU, SDAUN) is known. 


JSUBBR 

This is an abbreviated version of the surface cell 
list JSURF, including only those cells that are sub- 
divided. 

LTABBR (2, NSDELT) 

An abbreviated version of the element table LTABL, 
including only subdivided elts. An extra word is 
needed to indicate the location of elts. In 
LTABL (17,17,33) location is equivalent to position 
in the array. 
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SDYOU, SDAUN, SDPEE, SDARR, SDDIV, SDROUS 


These are the arrays used by TRILIN to .compute 
potentials at subdivided points . They correspond 
to the grid point arrays U, AUN, P, R, DIV, and 
ROUS, respectively. The ’’SD 1 ’ arrays are one-dimen- 
sional, as subdivided points are not arranged in an 
orderly three-dimensional fashion. 

PSURF (NSUB, NSUB , NUMFA) , PEP (NSUB, NUMED) 

These arrays are a more convenient way to reference 
the SDYOU array „ SDYOU begins with face subdivision 
points (PSURF) and ends with edge subdivision points. 
Consequently PSURF (1,1,1) is equivalenced to SDYOU (1) . 
PED (1,1) is equivalent to SDYOU (NSUB 2 • NUMFA + 1), 
the location of the first edge sd point. 

RSURFS (NSUB, NSUB, NUMFA) , REDS (NSUB, NUMED) 

Like PSURF and PED, but these arrays are equivalenced 
to locations in the SDAUN array. RSURFS (1,1,1) 
equiv SDAUN (1) ... REDS (1,1) equiv SDAUN (NSUB 2 

NUMFA + 1) . 
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7.12 SUBDIVISION SUBROUTINE SUMMARIES 


BAKINT 

Calling Sequence : SUBROUTINE BAKINT (TERP , NPT , BEG, XEND) 

Called By : EMPDI , EMPSLI 

Purpose : Linear back-interpolation. Values stored in TERP 

for one, two, or three equally spaced points are back- inter- 
polated to endpoints BEG and XEND. 

BILBAK 

Calling Sequence : SUBROUTINE BILBAK (RESARR, NSUB) 

Purpose : Bilinear back-interpolation. The values of the top 

face (Z =2) of a subdivided cube, stored in RESARR, are back- 
interpolated to the four corners. NSUB is the number of sub- 
divisions = 1, 2, or 3 . 

BILINT 

Calling Sequence : SUBROUTINE BILINT (VT, DIPOT, NSUB) 

Purpose : Bilinear interpolation. Vertex values from VT are 

interpolated for values on the top face of a subdivided cube. 
DIPOT stores the results. 

CLOKRO 

Calling Sequence : SUBROUTINE CLOKRO (VFROM, VTO) 

Purpose : Clockwise rotation of cube. Vertex values in VFROM 

are rotated 90° around an axis parallel to the Z axis, re- 
sulting in cube VTO. The rotation is clockwise looking from 
plus Z . 
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CORNAD 


Calling Sequence ; SUBROUTINE CORNAD (IJK, NORM, PIECE, 

BIGARR, NX, NY, NZ) 

Called By : GETDIV, EDGADD, CORNAD 

Purpose : Add to surface corners. Finds locations in BIGARR 

corresponding to corner nodes of square surface cell with 
origin IJK and normal NORM. Adds PIECE to each of these 
locations . 

COTELL 

Calling Sequence : SUBROUTINE COTELL (AUN, UP) 

Purpose : Prints out all non-zero values in two big arrays 

(17 x 17 x 17) . 

DIAGCA 

Calling Sequence : SUBROUTINE DIAGCA (DX, DY, DZ) 

Purpose : Diagonal calculation. Calculates diagonal term 

of coproduct for vertices of a given rectangular parallele- 
piped. It is 1/9 (DY-DZ/DX + DX*DZ/DY + DX*DY/DZ). 
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DICER 


Calling Sequence : SUBROUTINE DICER (VT, AUT, IABC, IJK, 

NSUB, DZ, NUMFA , LOFACE , PSURF, RSURFS , 
NUMED, LOEDGE , PED, REDS) 

Called By : SUBDIV, SURFSD 

Purpose ; Residual calculation for face-sd cells and elts. 
DICER retrieves U values, makes successive calls to RESID 
or TRIRES, and stores the resulting AUN values. 

Arguments include the vertex U values (VT) and the 
origin (IJK) and orientation (IABC) of the subdivided face. 

U values for subdivided points are found in PSURF and PED 
arrays . 

FILDI performs a bilinear interpolation for the top 

face. The subdivided object can now be treated as a set of 
2 

(HSUB + 1) rectangular parallelepipeds. 

TRIRES or RESID is called for each of the rectangular 
parallelepipeds . 

The residual (AUN values) on top are back-interpolated 
to the vertices by EMPDI. Residuals of sd nodes are added to 
RSURFS and REDS arrays. Vertex residuals are returned in 
AUT. 

(See also section "RESIDUALS AND SUBDIVISION".) 


EDGABC 

Calling Sequence : SUBROUTINE EDGABC (NSID, IFAABC, NSIABC) 

Purpose ; Generate orientation vector for an edge. NSIABC 
is returned as the orientation vector of side number NSID of 
the face with orientation vector IFAABC. 
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EDGADD 


Calling Sequence : SUBROUTINE EDGADD (IJK, NORM, PIECE, 

BIGARR, NX, NY, NZ , SMLARR, NUMFA, LOFACE , 
NUMED, LOEDGE) 

Called By : Q UP DAT , GETDIV 

Purpose ; Add to nodes of edge-sd surface. Locations cor- 
responding to sd points are found in SMLARR. Corner points 
are in BIGARR. Add in the following proportions: 

Edge-sd nodes: 2 * PIECE. 

Corner nodes on sd edge: 1 • PIECE. 

Corner nodes opposite sd edge: (NSUB + 1) • PIECE 

This distribution corresponds to uniformly distributed charge 
or back-interpolation from (NSUB + 1) equal rectangular 
parallelepipeds. 

EDGETE 

Calling Sequence : SUBROUTINE EDGETE (IXELT, IXEDGE, LTABLE , 

NX, NY, NZ, NSDELT , LTABBR) 

Called By : EDLTAB 

Purpose : Edge test for possible sd elt. Checks elt with 

origin IXELT to see if edge-sd. Other possibilities are 
that it may be filled or already sd. If already sd, a con- 
sistency check is made. If not already sd and empty, the 
element table is altered. 
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EDGIND 


Calling Sequence ; SUBROUTINE EDGIND (NSID, IABC, IJK, IS,Q) 

Purpose : Edge index code. For edge number NSID of the given 

face (IABC, IJK) , the edge index code (ISQ) is returned. 

EDLTAB 

calling Sequence : SUBROUTINE EDLTAB (IXFAC, LTABLE , NX, NY, 

NZ, NSDELT , LTABBR) 

Purpose : Edge-sd elt test. Given a sd face, test the eight 

elts (call EDGETE) that share one edge with this face. See 
if they are edge-sd. 

EMPDI 

Calling Sequence : SUBROUTINE EMPDI (RESARR, AUT, FACRES , 

EDGRES, NSUB) 

Called By : DICER 

Purpose : Break RESARR array into component parts. Calls 

BILBAK for back-interpolation of top (Z = 2) face of sub- 
divided cube. Puts vertex values into AUT, sd face into 
FACRES, and sd edge into EDGRES. 

Complements FILDI. 


EMPSLI 

Calling Sequence : SUBROUTINE EMPSLI (RESARR, AUT, EDGRES, 

NSUB) 

Called By : SLICE 

Purpose : Breaks single edge RESARR into vertex values (AUT) 

and one edge values (EDGRES) . 

Complements FILSLI. 
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EXPAND 


Calling Sequence : SUBROUTINE EXPAND (IX , IDIR, IJK) 

Purpose : Expand single word code into direction (IDIR = 

1, 2 or 3) and origin (IJK) . Inverse operation of SQUISH. 

FACADD 

Calling Sequence : SUBROUTINE FACADD (IJK, NORM, PIECE, 

BIGARR, NX, NY, NZ , SMLARR, NUMFA, LOFACE, 
NUMED, LOEDGE) 

Called By : QUPDAT , GETDIV 

Purpose : Add to nodes of face-sd surface. Locations of sd 

face and edge points in SMLARR, grid points in BIGARR. For 
even distribution, sd face nodes get 4 • PIECE, edge nodes 
get 2 * PIECE, corner nodes get 1 • PIECE. 

Similar to EDGADD. 


FACETE 

Calling Sequence : SUBROUTINE FACETE (LTELT , IJK, IABC, 

LTABLE , NX, NY, NZ , NSDELT , LTABBR) 

Called By : FALTAB 

Purpose : Test for face-sd elts. Alter element table if 

elt is empty and not .already sd. Otherwise test for con- 
sistency. 

Similar to EDGETE. 


FACIND 

Calling Sequence : SUBROUTINE FACIND (IABC, IJK, ISQ) 

Purpose : Index code for face. Construct one-word index code 

(ISQ) for face with orientation IABC on elt with origin IJK. 
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FALTAB 


Calling Sequence : SUBROUTINE FALTAB (IXORG, NORM, LTABLE, 

NX, NY, NZ, NSDELT, LTABBR) 

Purpose ; Face-sd elt test. Given an sd face, test (call 
FACETE) the elt in front and behind. One should be filled, 
the other face-sd. 

FDOT 

Calling Sequence ; FUNCTION FDOT (VONE, VTWO, LENGTH) 

P-urpose : Floating dot product. Take dot product of two real 

vectors of given length. 


FILDI 

Calling Sequence : SUBROUTINE FILDI (VT, FRECK, XNOTCH, DIPOT, 

NSUB) 

Called By : DICER 

Purpose ; Fill DIPOT (5,5,2) array from its component- parts — 
vertices (VT) , face (FRECK) , and edges (XNOTCH) . Top face is 
interpolated. Preparation for calls to residual routine. 

Complements EMPDI . 


FILLAR 

Calling Sequence : SUBROUTINE FILLAR (ICORN, P, FILLIN, PSURF, 

PED, LOFACE, LOEDGE) 

Purpose : Fill an array (FILLIN) with potential values for an 

area that includes sd node’s . 
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FILSLI 


Calling Sequence : SUBROUTINE FILSLI. (VT, XNOTCH, SLIPOT, 

NSUB) 

Called By : SLICE 

Purpose : Build SLIPOT (5,2,2) from vertices (VT) , one sub- 

divided edge (XNOTCH) , and three edge interpolation. 

Similar to FILDI . 

Complements EMPSLI. 


FINDED 

Calling Sequence : SUBROUTINE FINDED (IJK, NORM, NUMFA, 

LOFACE, NUMED, LOEDGE, IFOND , JSID) 

Purpose : Find an edge. Get storage location of the edge 

for a surface that is edge-sd. Will not work for face-sd 
cell . 

GENVER 

Calling Sequence : SUBROUTINE GENVER (IJK, NORM, I HAVER) 

Purpose : Generate vertices. Returns in IHAVER (3,4) the 

coordinates of the four corners of a surface with origin 
IJK and normal NORM. 

IDOT 

Calling Sequence : FUNCTION IDOT (JV, KV, NLEN) 

Purpose : Integer dot product. Value is dot product of two 

integer arrays with length NLEN. 
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INDEX 


Calling Sequence : FUNCTION INDEX (IDENT, IARRAY, LENGTH) 

Purpose ; Look up index number. IARRAY is an array of index 
codes of given length. If IDENT corresponds to one of the 
entries INDEX returns the entry number. If it is not on the 
list, value of INDEX is 0. 

LININT 

Calling Sequence ; SUBROUTINE LININT (BEG, XEND, NPT , TERP) 

Purpose ; Linear interpolation BEG and XEND are the values at 
the endpoints . TERP returns with values for NPT equally 
spaced intermediate points (NPT <_ 3) . 

LONEDG 

Calling Sequence ; SUBROUTINE LONEDG (IJK, IOR, IXED) 

Purpose ; Edge index code for one sd edge. Given elt origin 
(IJK) and sd edge orientation (IOR) returns edge index code 
IXED. 

MAFILL 

Calling Sequence ; FUNCTION MAFILL (LTAB) 

Purpose ; Filled elt test. MAFILL ^ 0 if element table entry 
LTAB refers to a completely filled element. 

MATMUL 

Calling Sequence ; SUBROUTINE MATMUL (ARRAY, VECTOR, RESULT, 

NDIM) 

Purpose ; Matrix multiplication. The square array is multi- 
plied by the vector to yield the result. All are real-valued. 
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OFFSET 


Calling Sequence : SUBROUTINE OFFSET (IABC, IOFF) 

purpose ; Compute offset (IOFF) from elt origin corresponding 
to orientation vector IABC. 

PTAA, PTAB, PTAC 

Calling Sequence : SUBROUTINE PTAC (LABEL, ID, JD, KD, ARRAY, 

IR, JR, KR) 

Purpose : Print array. PTAA for one-dimensional arrays, 

PTAB for two-dimensional and PTAC for three-dimensional. 

The array is dimensioned (ID, JD, KD) . Values are printed 
for X = 1, IR . . . Y = 1, JR ... Z =1, KR. The label is 
printed first. 

PTAQ 

Calling Sequence : 

Purpose : Print a 

sent the vertices 
format suggesting 

RESID 

Calling Sequence : SUBROUTINE RESID (SKIN, SKOUT, AL, S) 

Called By : DICER, SLICE 

Purpose : Calculate residuals for subsection of a volume elt. 

SKIN holds U values for the vertices of a rectangular parallele- 
piped with dx = dy = S, dz = AL. SKOUT returns corresponding 
AUN values. 

TRIRES performs the same function for surface cells. 


SUBROUTINE PTAQ (LABEL, Q) 

cube. The array Q (2,2,2) is taken to repre- 
of a cube. The values are printed in a 
a perspective drawing of a cube. 
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SDINPU 


Calling Sequence : SUBROUTINE SDINPU (no arguments) 

Called By ; INPUT 

Purpose : Input subdivision information. Takes user specifi- 

cation for subdivided area, and alters surface cell list 
(JSURF) for face-sd cells. Sets NSUB. 

SDLIST 

Calling Sequence : SUBROUTINE SDLIST (LTABLE, NX, NY, NZ , 

IOBJIN) 

Purpose : Alter lists and make new lists for subdivision. 

Surface cell list and element table altered for sd cells and 
elts . Forms lists of sd faces, sd edges, sd cells, and sd 
elts. Some of the work done through calls to FALTAB and 
EDLTAB. 

SLICE 

Calling Sequence : SUBROUTINE SLICE (UT, AUT, IABC, IJK, NSUB, 

DZ, NUMFA, LOFACE, PSURF, RSURFS, NUMED, 
LOEDGE, PED, REDS) 

Called By : SUBDIV, SURFSD 

Purpose : Residual calculation for edge-sd cells and elts. 

Structure is parallel to DICER (see DICER) . Difference is 
that only one edge has sd nodes and three edges must be 
interpolated (see section "RESIDUALS AND SUBDIVISION"). 
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SPILSD 


Calling Sequence : SUBROUTINE SPILSD (ICORN, IDELTS, P, 

FILLIN, PSURF, PED, LOFACE , LOEDGE) 

Purpose ; Output. Prints potential values for given 
rectangular area including some subdivided cell(s). 


SQUISH 

Calling Sequence ; SUBROUTINE SQUISH (IDIR, IJK, ISQ) 

Purpose ; Produce index code ISQ. The four numbers IDIR 
and IJK (3) are encoded into one word. 

Inverse operation to EXPAND. 


SUBDIV 

Calling Sequence ; SUBROUTINE SUBDIV (UP , AUN, NX, NY, NZ , 

DOTPRO, SDYOU, SDAUN , NUMFA, LOFACE, 
NUMED, LOEDGE, NSDELT , LTABBR) 

Called By ; ISPACE 

Purpose ; Subdivided elt residual calculation, top level 
routine. Cycle over sd elts. For each retrieve vertex U 
values and transform to standard orientation. Call DICER 
or SLICE. Add resulting vertex AUN values (residuals) to 
grid node AUN array. 

Similar to SURFSD. 
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SUBSGR 


Calling Sequence : SUBROUTINE SUBSCR (ICORN, LH, MV, MYORG, 

LITLH, LITMV) 

Called By ; FILLAR 

Purpose ; Subscript computation for data organizing routine 
FILLAR. Given offsets of a single square cell in area of 
interest, returns grid coordinates of square origin and 
FILLAR coordinates for sd points. 


SURFSD 

Calling Sequence ; SUBROUTINE SURFSD (UP, AUN, NX, NY, NZ , 

UPCOND, AUCON, SDYOU, SDAUN, NUMFA, LOFACE, 
NUMED, LOEDGE, NSDURF, JSUBBR) 

Called By ; ISPACE 

Purpose ; Sd surface cell residual calculation, top level 
routine. Cycles over sd surface cells. For each, retrieves 
U values for four corners and underlying conductor. Trans- 
lates edge-sd cell into standard orientation. Calls DICER 
or SLICE. Adds grid node AUN values (residuals) to AUN 
array. 

Similar to SUBDIV. 


TRANS D 

Calling Sequence ; SUBROUTINE TRANSD (VFROM, VTO, NS) 
Called By ; SURFSD 

Purpose ; Transform edge-sd cells. Sd edge is on bottom 
(minus Z) face. It is side number 1, 2, 3, or 4. TRANSD 
transforms the vertices to bring sd edge to side 1. Also 
performs inverse transformation. 
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TRIRES 


Calling Sequence : SUBROUTINE TRIRES (SKIN,SKOUT, DX, DY, DZ) 

Called By ; DICER, SLICE 

Purpose : Residual calculation for rectangular parallelepiped 

section of surface cell. DZ is the effective thickness of 
surface cell. DX and DY depend on NSUB and whether surface 
cell is face-sd or edge-sd. Performs appropriate matrix 
multiplications with a damping term' and scalar multiplication 
by the dielectric constant for this cell. 

Similar to RESID which handles volume elts. 

UNSURF 

Calling Sequence : SUBROUTINE UNSURF (JSUWOR, MCONDU, IJK, 

NORM, MATL, DUMA, DUMB) 

Purpose : Decodes entry in surface cell list (JSURF) for in~ 

formation about a single cell. 
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7 . 13 CODE STRUCTURE 


The following diagrams illustrate the new subdivision 
routines and old routines that have been altered for sub- 
division. For the sake of clarity, some routines are included 
even though they have not been changed. These are marked with 
an asterisk (*) . 
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LEVEL 1 



Not changed by subdivision. 
































LEVEL 3 




Not changed by subdivision. 
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8 . PARTICLE DETECTORS 


NASCAP contains a set of routines which permit the 
user to obtain a fine resolution calculation of the energy 
flux density function as observed at user-specified locations 
on the satellite. The routines have been constructed so that 
the calculations render information in a form which resembles 
that which might be obtained experimentally from a particle 
detector located on the satellite. In order to use the rou- 
tines the NASCAP user must specify that a particle detector 
is located upon one or more surface cells of the satellite 
model. In terms of a spherical coordinate system which is 
located upon the surface cell the user may then elect to have 
the energy flux density (as observed by the detector) plotted 
as a function of any one of the three variables: detector 

polar angle, detector azimuthal angle, and kinetic energy of 
the incident particles. One variable is selected as the 
independent variable to be swept through a user-specified 
range of values while the other two remain fixed at user- 
selected values. 

In addition, the user may provide a description of the 
particle detector's resolution. The calculated value of the 
energy flux density is then obtained as an integral average 
over the detector's angular and energy apertures. 

The energy flux density for each detector which the 
user has "activated" is calculated and plotted as a function 
of the selected independent variable at user-selected points 
within the NASCAP runstream. Plots for both electron and 
proton energy flux density are combined on a single overlay 
plot. A separate plot is made for each detector. 

Also available to the user is an option to plot the 
particle trajectories generated during reverse-trajectory 
evaluation of the energy flux density function. If this 
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option is selected then each time a detector flux plot is 
generated, a set of particle plots is also generated. In 
order to constrain the number of plotted trajectories to a 
reasonable number, only one particle trajectory is plotted 
for each value of the selected independent variable at which 
the detector energy flux density integral is evaluated even 
though more particle trajectories may have been generated. 
(Referring to Eq. (8.10) in the section entitled "Detector 
Energy Flux Density Measurement", the trajectory which gets 
plotted is the one for which i = j = k - 1.) Separate plots 
of particle trajectories projected onto the Y-Z, X-Z and 
X-Y planes are produced for each of the two particle species 
electrons and protons. Thus six particle plots are produced 
for each detector at each multiple of the NASCAP time cycle 
step selected. (Note the potential for an excessive number 
of plots if care is not exercised!) Each particle plot also 
includes a silhouette of the satellite. 
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8.1 


SPECIFICATION OF DETECTOR ORIENTATION 


In order to specify the direction in which a particle 
detector is to measure energy flux density it is necessary 
for the NASCAP user to be familiar with three different right- 
handed coordinate systems. 

1 . The Satellite Coordinate System 

This is the coordinate system in which the satellite 

building blocks are defined. 

2 . The Surface Cell Coordinate System 

This system is obtained by performing the following 

operations on the satellite coordinate system : 

a. The X-Y plane is rotated counterclockwise through 
an angle ® about the Z-axis until the +X-axis co- 
incides with the X-Y -plane projection of the sur- 
face cell normal vector. 

b. The Z-X plane is rotated counterclockwise through 
an angle if; about the Y-axis until the +Z-axis co- 
incides with the cell normal vector. 

c. The origin of the system (obtained by successive 
rotations in a. and b. is translated to the center 
of the surface cell face. 

The procedure of a - b is illustrated in the following 
diagram. The double-primed coordinate axes are for the cell 
system while the unprimed coordinate axes are for the satel- 
lite system. Table 8.1 lists \p and ® for each of ‘the 26 pos- 
sible cell normals which NASCAP can produce. 
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Table 8. 

1 . 

Cell Coordinate System 
Satellite System 

Rotation With Respect 

Cell 

Normal 

Type 

i 

$ 

-1 

-1 

-1 


125.26° 

225.00° 

-1 

-I 

0 


90.00 

225.00 

-1 

-1 

1 


54.74 

225.00 

-1 

0 

-1 


135.00 

180.00 

-1 

0 

0 


90.00 

180.00 

-1 

0 

1 


45.00- 

180.00 

-1 

1 

-1 


125.26 

135.00 

-1 

1 

0 


90.00 

135.00 

-1 

1 

1 


54.74 

135.00 

0 

-1 

-1 


135.00 

270.00 

0 

-1 

0 


90.00 

270.00 

0 

-I 

1 


45.00 

270.00 

0 

0 

-1 


180.00 

0.00 

0 

0 

0 


- 

- 

0 

0 

1 


0.00 

0.00 

0 

1 

-1 


135.00 

90.00 

0 

1 

0 


90.00 

90.00 

0 

1 

1 


45.00 

90.00 

1 

-1 

-1 


125.26 

315.00 

1 

-1 

0 


90.00 

315.00 

1 

-1 

1 


54.74 

315.00 

' 1 

0 

-1 


135.00 

0.00 

1 

0 

0 


90.00 

0.00 

1 

0 

1 


45.00 

0.00 

1 

1 

-1 


125.26 

45.00 

1 

. 1 

0 


90.00 

45.00 

1 

1 

1 


54.74 

45.00 
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The surface cell coordinate system is the coordinate 
system in which the NASCAP user can specify the direction 
in which the detector is to be pointed. This is done most 
conveniently by transforming the rectangular surface cell 
coordinate system into a spherical one. The detector orien- 
tation is then specified by <J> , the azimuthal angle and 0, 
the polar angle. The following table contains a few key 
combinations that may help to clarify the use of <f> and 0 . 


4 > 


0 


Direction Detector Points 


0° 0° +2 axis of cell system (cell normal) 


0° 

90° 

+X 

axis 

of 

cell 

system 

90° 

90° 

+Y 

axis 

of 

cell 

system 

180° 

90° 

-X 

axis 

of 

cell 

system 

270° 

90° 

-Y 

axis 

of 

cell 

system 


When measuring the energy flux density function it is 
probably the energy dependence of the function which will 
be of interest so that 9 and <j> will be maintained at 'fixed 
values. It 'is possible, however, to display the energy flux 
density as a function of either 0 or <J> for a user-specified 
angular range. This might be a desirable thing to do in 
cases for which the detector does not have a clear view in 
all directions (i.e., the detector might be shadowed in some 
directions by other parts of the satellite) . 
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8.2 


THE DETECTOR COORDINATE SYSTEM 


This system is obtained by performing the following 
operations on the surface cell coordinate system. Let v be 
the vector which points in the direction in which the detector 
energy flux density is to be measured. (v has coordinates 
<j> , 6 in the spherical coordinate system of the surface cell.) 

a. The X-Y plane is rotated counterclockwise through 
the angle <p about the Z axis until the +X axis 
coincides with the projection of v onto the X-Y 
plane of the surface cell coordinate system. 

b. The Z-X plane is rotated counterclockwise through 

the angle 0 about the Y axis until the +Z axis 

• « « ^ 
coincides with v. 

The NASCAP user will normally not be directly concerned 
with the details of how the detector coordinate system relates 
to the other two coordinate systems . All that one need be 
aware of is that the +Z axis of the detector coordinate system 
is the central axis of the hemispherical cap which is used to 
model the angular aperture of the detector. (See Section 8.4, 
"Detector Energy Flux Density Measurement".) 
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8.3 


CALCULATION OF ENERGY FLUX AT A CELL SURFACE 


In order to obtain an expression for the energy flux 
density measured by a detector located at a given surface 
cell of the satellite model it is helpful to first consider 
the general problem of calculating the total energy flux 
which is incident at the surface of the cell due to the am- 
bient plasma environment. Let k be the unit normal vector 
for the surface cell. Using the cell’s rectangular coordinate 
system (obtained by appropriate rotation of the satellite co- 

A 

ordinate system) with the +Z axis in the direction of k, the 
energy flux at the cell surface center is calculated as fol- 
lows : 


03 CO 00 / > | 2 

m i V ol ■ 


«• ■ - £ / L J 


2e 


1 

e 


(eV -k) f (V J d V 


o o 


— CO —CO —CO 


■± z 

where t = energy flux vector (eV/(M -sec)) at cell surface 


^ = vi + vj+vk (velocity at surface 

o x y z 

d~V = dv dv dv 
o x y z 

—y 

f (V ) = phase space density function evaluated at the 
surface for particles with velocity V . 

m = mass of incident particles. 


e = charge of incident particles (electrons or 
protons) . 

To evaluate the integral it is expedient to change 
from rectangular to spherical coordinates. The necessary 
substitutions are: 
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d 2 V = v 2 sin0 d6 dd> dv 
o o r o 

-V A 

(eV *k) = ev cos0 
o o 

With- these the energy flux integral transforms into 








sin0 d6 difr dv 
o r o 


« 2ir 1 





f 0 (w,4.,v 0 ) 



d<J> dv 

o 


where u = cos© . 

Next we change variables from v q {the magnitude of 
the velocity at the surface) to E q (the kinetic energy at 
the surface of the cell) . 

Let 
1 2 

7 T m v = e E (factor e because E is in eV) 

2 o o' o 

edE 

dv = 

o mv 

o 

Then the energy flux integral becomes 


“ 2ir 1 


o “ - £ / / / 2 (s) E o 


f 0 (U/^,E 0 ) _ dd, dE o 


0 "0 "0 


00 2 tt 

-s// 


e \ 2 2 

2 - E f (E ,0 ) u dfi dE 
’ m / o o o o o 


0 0 


( 8 . 1 ) 
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where we have introduced the solid angle fi as an integration 
variable for notational convenience. 

Finally we introduce the energy flux density function 
G(E,fi). From the definition we know that 






G (E ) V dfi dE 
o o' o o 


( 8 . 2 ) 


Comparison of the integrands in Eqs. (8.1) and (8.2) 
yields the identity 



f (E 
o o 



G o (E o' D o> 


(8.3) 


This equation provides the key to correct evaluation 
of G at the cell surface by reverse trajectory particle 
tracking. Since f is a constant along a particle trajectory 
we have for particles emitted at the cell surface with initial 

♦ • • i 

energy E q and velocity vector m the direction of 


f (E . ^ ) = f (E ,$ ) 

O O O CO OD ' CO 


(8.4) 


G o (E o'°o> 


E‘ 


G (E ) 

CO 00 ' rrt ' 


E 


(8.5) 


where 


f = phase space density function 
surface. 

f m = phase space density function 

G q = energy flux density function 
surface. 


evaluated at cell 

evaluated at infinity, 
evaluated at cell 
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= energy flux density function evaluated at infinity. 

E Q - initial kinetic energy of particle to be tracked. 

= kinetic energy of particle after reverse trajectory 
tracking to infinity. 

= initial velocity direction vector of particle at 
cell surface. 

^co “ velocity direction vector for particle after • 
reverse trajectory tracking to infinity. 

Therefore, if the energy flux density function G is 
known at infinity then using reverse trajectory particle 
tracking the energy flux at the surface of the cell may be 
computed from 



This equation is used (with slight modification) in 
the following section to arrive at an expression for the 
energy flux density which is measured by a particle detector 
located on the surface cell. 
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8.4 


DETECTOR ENERGY FLUX DENSITY MEASUREMENT 


Consider an ideal particle detector located at the 

center of a surface cell and oriented such that the +2 axis 

of the detector* s rectangular coordinate system points in the 

direction in which the energy flux density is to be measured. 

Assume that the detector’s rectangular coordinate system is 

■)* , 

transformed into a spherical coordinate system and that n is 
the unit normal pointing in the direction in which the 
particle detector is pointing. If the detector responds to 
particles having energy E q then the measured value for the 
energy flux density will be 


BE 

o 


- s e o{ 2 (s-r f o< E o' s >> 


-n E 2 
o 


G (E ) 

CO ' CO ' CO ' 


( 8 . 7 ) 


where 

B q - energy flux vector at cell surface. 

f Q = phase space density function evaluated at cell 
surface. 

E q = kinetic energy (in eV) of particles at cell sur- 
face. 

<f> = azimuthal angle in detector coordinate system. 

0 — polar angle in detector coordinate system. 

E^ = kinetic energy (in eV) of particles after reverse 
trajectory tracking to infinity. 

= velocity direction vector for particles after 
reverse trajectory tracking to infinity. 

m = mass of particles. 
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G^CE,.^) = energy flux density function (i.e. , tabulated 

Deforest data, for example. Units are eV/ 

2 

(sec-cm -sr-eV).) 


In practice, a real detector responds to a finite range 

of particle energies and velocity vectors so that the energy 

flux density observed by the detector is really 9 2 *f /9E 9fi 

o o 

averaged over the energy and angular apertures of the detector. 
Thus, in general, the value which the detector yields is of 
the form 


9 2 ~E 


9E 9ft 
o 



n 2 i co) 

W < E o' a> E o —2 


cosQ dft dE 

o 



W<E o ,«) dft dE Q 


( 8 . 8 ) 


where W(E Q ,fi) is a weight ' function which describes the charac- 
teristics of the detector’s energy and angular apertures. 

In the current version of NASCAP it has been assumed 
that the energy aperture of a detector is a rectangular 
weight function which has a value of 1 from E to E + AE and 
0 elsewhere. The angular aperture is assumed to be a hemi- 
spherical cap of "width" A0 in the polar angle. {The vector 
n passes through the center of this cap.) The weight has a 
value of 1 anywhere on the cap and 0 elsewhere. Therefore, 
NASCAP detectors compute the energy flux density from 
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3E 3ft 
o 


o _ E 


E+AE 2tt 1 

" S / / / 

*Tn «/q •'nr 


E 


2 

cos A0 


(G (E ) ) 

I CO ' CO 9 CO 7 I 

~ £ J 


T- a f dE t 


E+AE 2tt A0 

/ // 

^•n •'n • / n 


sin0 d0 <3.(|) dE 


E+AE 

f 

2ir 1 i 

r r b 2 


[ du 2 

l J 

1 J 2 ° I 

0 •'cos A0 

1 E 2 1 

l. CO J 

1 2 


dcj) dE 


(1 - cosA0} 2 tt AE 


(3.9) 


where \x = cos0. 

This integral is evaluated by a three-dimensional ap- 
proximation formula which uses the mid-point rule with n^ 

points to compute the integral over <|>,-the mid-point rule 

2 

with points to compute the integral over u , and the 
even-order Gauss-Legendre formula with n points to compute 
the integral over E q . The composite formula used^is: 


) 2 E 


n. 


3E 3 n 
o 


“~n< 


$ 2 
2u \ V- (1-cos A0 ) 

n . 2 -j 2n 


n 




n /2 
z , 


♦' “i 


<i> 


E(“)E |v? F 

j=l 




k=l 


+ W. k e! k J 


1 


(-l-cosA0) 2tt Ae 


n, n. 

<f> $ 


n /2 

e , 


= -ny 


L E EM +w - k 4 

i=l j=l k=l ^ 


(8.10) 
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where 


e k ■ “ <X* - 1) + E 
♦i ■ 57 (i - 1/2) 

lij = Ay 2 (j - 1/2) + (1 - n Ay 2 ) denotes y 2 } 

and 

Ay 2 = (1 - cos 2 A8)/n^ 


(1 + cosA0) , ... , , . _-4 . , . , -2 , 

y = - — - (multiply by 10 no put units in cm ) 

4 n , n . 

\p q> 


F( *i' u j' e k ) = 


G (E ,S5 ) 

CO CO 1 CO 


E 


(returned by FSPACE) 


ijk 


E^ and are the final energy and velocity vector 

respectively of a particle after reverse trajectory tracking 

from the center of the surface cell (beginning .with initial 

2 

velocity specified by © . , y., and £,)• to infinity. 

1 0 K 

The and are the Gauss-Legendre integration co- 
efficients for n an even integer. (Note that x^ = an< ^ 

= W_^,.) A slightly modified formula is used to permit 
n = 1. (Also note that -1 < Xv < 1 for all k.) 

It should be noted that although the detector energy 
flux integral includes only contributions from the ambient 
plasma environment it is possible that some particle tra- 
jectories will yield < 0. This could occur if the particle 
originates from another part of the satellite, for example. 
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For detector particle trajectory plotting purposes all parti- 
cles must be tracked regardless of origin. Therefore, a test 
is made within the innermost integral summation loop to 
determine if E >0. If it is not then no attempt is made to 

CO — — 

evaluate G and G is assumed to be 0. 
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8.5 


NASCAP PARTICLE DETECTOR ACTIVATION 


Each of up to twenty particle detectors is activated by 
the appearance of a detector keyword definition sequence in 
the NASCAP runstream. A detector keyword definition sequence 
consists of a "DETECT" keyword card which is optionally followed 
by one or more of the keyword cards to be described below. A 
detector keyword definition sequence is terminated either by 
an "END" keyword card or else by another detector keyword 
definition sequence. If more than one detector is to be acti- 
vated "simultaneously" the associated keyword definition 
sequences must be consecutive. The last detector keyword 
definition sequence to appear on the NASCAP keyword file must 
always be terminated by an "END" card. In addition, if no 
other NASCAP runstream options are to follow the detect key- 
word sequences a second "END" card should follow the detect 
"END" card. The reason for this is that the first "END" 
terminates the detector run while the last "END" terminates 
the NASCAP runstream. 

\ 

If it is desired to call the detector routine more 
than once from the NASCAP runstream using the same set of de- 
tector options then the detector keyword definition sequence (s) 
should reside upon its (their) own separate file. 

If the FORTRAN file number used is 23, for example, 
then each time a DETECT 23 card appears in the NASCAP run- 
stream detector routine execution will be initiated using 
the detector keyword options in file 23. If the detector key- 
word options reside in a file other than the NASCAP runstream 
then the DETECT card on the runstream should only be followed 
by an "END" card if no other functions are to be performed 
by NASCAP following detection. As an example, assume the 
detector keyword file is 23 and that following detection, we 
desire to do HIDCEL followed by TRILIN and finally another 
detect. The NASCAP runstream would look like: 
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DETECT 23 
HIDCEL 
TRILIN 
DETECT 23 
END 

If, on the other hand, the detector keywords are to be inserted 
in the N&SCAP runstream itself then the runstream might look 
something like: 



END 


first detector run 


second detector run 
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8.6 


NASCAP DETECTOR SCRATCH FILES 


The NASCAP detector routines require the use of two 
scratch files. The NASCAP file IAUN is used as a temporary 
file for storing detector flux measurement information. In 
general, the user need not concern himself with this. The 
other file used is the NASCAP file IPART.- It is the user's 
responsibility to make sure this file has been assigned for 
scratch usage. Since the file is used for particle trajectory 
plotting it is recommended that additional space beyond the 
computing system default for normal files be allocated. The 
default value of IPART is -28. (The negative value is re- 
quired to indicate that an old NASCAP routine, PARPLT, is 
not to be activated. The absolute value of IPART is used by 
the detector routines, however.) On the UNIVAC 1100/81 machine 
it is recommended that the runstream include the card 

@ASG , T 28 ,F///1000 
prior to the @XQT card. 
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8.7 


DETECTOR KEYWORD INPUT FILE 


At the beginning of each detector keyword sequence 
(appearance of the "DETECT" card) all parameters which describe 
the properties of the detector being activated assume a set 
of default values. The NASCAP user then has the option of 
changing any of these values to suit his requirements. This 
is accomplished by including one or more of the cards to be 
described below. The contents of each card consists of a 
mnemonic keyword, left justified in card columns 1-6, and 
possibly the value of a data variable associated with the 
keyword. The type of data required (if any) is determined by 
applying the standard FORTRAN convention for variable types 
to the keyword. Variables must appear in columns 21-30 of the 
keyword card. The format used is 110 for INTEGER type and 
F10.0 for REAL type. In a few cases an option may be speci- 
fied by a card having the form keyword = option where option 
is another mnemonic name. This type of card requires that 
there be no embedded blanks between the keyword, equal sign, 
and option. 
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8.8 


DESCRIPTION OF DETECTOR PARAMETER OPTIONS BY KEYWORD 


8 . 8.1 

ICELL 

COMMEN 

8 . 8.2 

NP 

NMU2 

NE 

NSTP 

DTH 


General Detector Definition Parameters 


Index number of the surface cell at which a detector 
is to be, activated. The detector will be placed at 
the center of this surface cell. The acceptable range 
is 1 £ ICELL £ 1250. The default value is ICELL = 1. 
(It should be remarked that the same value of ICELL 
may be specified in more than one detector keyword 
sequence so that it is possible to simulate several 
types of detectors, each located on the same surface 
all within the same run.) 

This card is for programmer convenience in identify- 
ing the purpose of keyword cards in the detector 
keyword input file. If this card appears, columns 
7-72 may contain any comment the user desires. It 
is echoed back on the keyword listing at execution 
time but is otherwise ignored. 


Detector Aperture Definition Parameters 


Number of points used for integration over the 
detector aperture azimuthal angle (n^ in Eq. (8.10). 
Acceptable range is NP >_ 1. Default is NP = 1. 

Number of points used for integration over the 
transformed detector aperture polar angle para- 
meter ( n j_n Eq. (8.10) . Acceptable range is 
NMU2 >_ 1. Default is NMU2 = 1. 

Number of points used for integration over the de- 
tector energy aperture (n £ in Eq. (8.10). Acceptable 
range is NE = 1 or 2 <_ NE <_ 12 where NE must be an 
even integer. Default is NE = 1. 

Maximum number of steps per particle allowed during 
reverse trajectory tracking. Acceptable range is 
1 < NSTP < 500. Default is NSTP = 100. In 
practice, values around 300 or more are likely to 
be necessary. 

Polar angular width of the hemispherical cap used 
to model the detector angular aperture (A6 in Eq. 
(8.10)'. Acceptable range is 0.0° £ DTH £ 90.0°. 

The default is DTH = 0.0°. 
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DE 

8.8.3 

INDVAR 


N 


ENERGY 

PHI 

THETA 


Width of detector energy aperture in eV (AE in 
Eq. (.8.10)). Acceptable range is DE >_ 0.0. De- 
fault is DE = 0.0." 


Independent Variables and Fixed Parameters- 


Keyword = option type variable used to select inde- 
pendent variable for detector energy flux density 
plots . The acceptable options and their results are 


Option 

ENERGY 

PHI 

THETA 


Independent 

Variable 

E 

<f> 

0 


Default is INDVAR=ENERGY . 


Independent Variable 
Axis Scale Type 

Log 

Linear 

Linear 


Number of increments used for the independent 
variable (number of points on the detector flux 
plot horizontal-axis.) Acceptable range is 
3 < N <_ 500. Default is N = 100. 

Any one of the following three keywords (para- 
meters) may be selected as the independent varia- 
ble for detector flux plots (see INDVAR above) . If 
a parameter is selected as the independent varia- 
ble then the value specified by the associated 
keyword below is the starting value of the parameter 
(The ending value is selected using the keyword 
FINALV.) Otherwise the parameter is held fixed at 
the value specified by the associated keyword card 
for all points which the independent variable 
assumes . 


Kinetic energy (in eV) of particles incident at 
the detector (E in Eq. (8.10)). Acceptable range is 
ENERGY > 0.0. Default value is ENERGY = 10.0. 

Azimuthal angle <f> of the detector (in degrees) as 
measured in the spherical coordinate system o£ the 
cell at which the detector is located. Acceptable 
range is 0.0° <_ PHI <_ +360.0°. Default is PHI = 0.0. 

Polar angle 0 of the detector (in degrees) as 
measured in the spherical coordinate system of the 
cell at which the detector is located. Acceptable 
range is -9 0.0° <_ THETA < +90.0°. Default is 
THETA = 0.0. 
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FINALV 


8.8.4 ' 
PSCALE 

LWPEN 

FLXMIN 

FLXMAX 

AUTOS 


(Note: Letting PHI, THETA, and DTH assume their de- 

fault values will produce a detector which looks along 
the normal to the surface cell ICELL upon which it is 
located. ) 

Final value (in appropriate units) of the independent 
variable selected by INDVAR. If INDVAR=THETA then 
acceptable range is -90.0° <_ FINALV <_ +90. .0°. If 
INDVAR=PHI acceptable range is FINALV <_ 720°. If 
IND VAR=ENE RG Y acceptable value is ENERGY < FINALV < 
50,000 eV. Default value is 4.999 x 10 4 eV (assumes 
default INDVAR=ENERGY was used) . If INDVAR ^ ENERGY 
then the user must explicitly input a value for FINALV . 


Plot Scaling Options 


Proton to electron energy flux density scale factor. 
This factor determines the separation of the proton 
and electron flux density curves on the overlay plot 
generated by the detector routines. The proton flux 
values are multiplied by the factor 10.0**PSCALE be- 
fore plotting. The acceptable range is PSCALE > 0.0. 
Default is PSCALE = 5. 

Line width of pen used to draw the proton flux curve. 
The default is LWPEN = 3 raster increments. (The 
electron flux curve is always drawn with a line 
1 raster increment wide.) The acceptable range is 
1 < LWPEN < 10 . 

Minimum value for the logarithmic energy flux den- 
sity scale on detector plots. (The scale is in 
units of eV/ (cm2-sec-sr-eV) . ) The acceptable range 
is FLXMIN > 0.0. Default is FLXMIN = 10 4 . (This 
parameter is ignored if the AUTOS option is in ef- 
fect. ) 

Maximum value for the logarithmic energy flux den- 
sity scale on detector plots. The acceptable range 
is FLXMAX > FLXMIN. Default is FLXMAX = 10 12 . (This 
parameter is ignored if the AUTOS option is in ef- 
fect . ) 

If this card is included the detector routines will 
automatically select the scale limits FLXMIN and 
FLXMAX at execution time so that a reasonable por- 
tion of the data is displayed. If this card appears 
then any values specified for FLXMIN and FLXMAX are 
ignored. The default is manual scaling by the user 
as specified by FLXMIN and FLXMAX. The AUTOS option 
is highly recommended, however. 
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8.8.5 Reverse Trajectory Particle Tracking 


PLPART 


VCODE 


8 . 8.6 


If this card is included in a detector keyword 
definition sequence then particle trajectory plots 
will be produced for the detector. 

Tolerance limit for the maximum distance in code 
units which particles are permitted to move at the 
first timestep after emission from the detector 
cell. Acceptable range is 0.0 < VCODE <_ 10'. 0 . The 
default is VCODE = 0.3 inner grid units . 


Specification of Environment for Detectors 

The detector routines require that the environment flux 
definition file be present. The only acceptable flux 
types are TYPE 3, Maxwell or DeForest. If NASCAP is 
using some other flux type than this then a special 
flux file must be prepared and activated by doing a 
RDOPT and setting IFLUX to the detector flux file unit 
number just prior to doing DETECT. Following DETECT 
it may be necessary to do RDOPT to reset IFLUX to the 
correct file number to be used for further computations 
performed by TRILIN . 
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9 . PARTICLE EMITTERS 


9.1 LOW DENSITY PARTICLE EMITTERS: GENERAL DESCRIPTION 

NASCAP has the capability to simulate satellite charging 
effects which result from one or more low density particle 
emitters placed at user-specified cell locations. Given that 
a particle emitter is to reside upon a particular surface cell 
of the satellite model the NASCAP user may, via keyword input, 
specify a number of parameters which describe the emitter 1 s 
characteristics. Included among these are: type of particles 

emitted (-electrons or protons) , total emission current, and 
direction of emission (specified in terms of the polar and 
azimuthal angles in a spherical coordinate system which is 
located upon the 'emitter 1 s surface cell ) . In addition, the 
user may choose one of several current density functions to 
represent the energy and angular characteristics of the emit- 
ter gun. If desired, the user may elect to have the trajec- 
tories which result from emitter gun forward particle tracking 
plotted at specified multiples of the NASCAP time cycle. 

(The satellite surface cell currents are, of course, corrected 
to include emitter current contributions at each NASCAP time- 
step whether or not plots are produced.) Separate plots of 
particle trajectories projected onto the X-Y, Y-Z , and Z-X 
planes are produced for each of up to five particle emitters 
which the user may have activated. Optionally, particle tra- 
jectories for all emitters may be combined on a single set 
consisting of three two-dimensional projection views. The user 
may also elect to have the three projection views plotted more 
than once, each time using a different maximum grid boundary. 
This permits the fine details of particles returning to the 
satellite to be viewed while also displaying the details of 
particles circling in large radii due. to the presence of a 
magnetic field, for example. Each particle trajectory plot 
also includes a silhouette projection of the satellite. 
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With a few minor exceptions, the sets of parameters 
which describe the characteristics of the emitters may be 
specified independently of one another. Thus one could, for 
example, define an electron emitter with a beam current of 
1 ya on surface cell 1 and a proton emitter with a beam cur- 
rent of 10 iaa on surface cell 470. The conductors underlying 
these two cells must be different, however. (See "Restric- 
tions " section . ) 
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9.2 


SPECIFICATION OF EMITTER BEAM ORIENTATION 


Particle emitter beam orientation is defined in vir- 
tually the same way as detector orientation. (See the section 
entitled "Specification of Detector Orientation".) One dif- 
ference to be noted' is that in the case of a detector, energy 
flux may be displayed as a function of 8 or <f> by "sweeping" 
one variable or the other through a specified range while the 
other remains fixed. It is assumed that an emitter does not 
change its orientation during satellite charging calculations, 
however, so that 9 and <b always remain fixed at user-selected 
values for each emitter. It should also be noted that the 
detector coordinate system (which shall now be referred to as 
the emitter coordinate system) is the coordinate system in 
which the integral of the emitter’s current density function 
is performed. 
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9.3 PARTICLE EMITTER BEAM CURRENT REPRESENTATION 

The emission current density functions of NASCAP parti- 
cle emitters contain a single energy spectral peak and are non 
zero only over a finite spatial range. In particular the emit 
ter current is envisioned as flowing out from the geometric 
center of the surface cell upon which the emitter is located 
and through a small hemispherical cap of width A0 in the polar 
angle. (The Z-axis of the emitter coordinate system passes 
through the center of this cap.) NASCAP generates a finite 
representation of the emission current by emitting a discrete 
set of particles. The angles and energies with which parti- 
cles are emitted are chosen in an optimized manner so that 

the total emission current is divided equally among each 

* 

particle of the set. Each particle is "pushed" in the 

"k i? 1 

electrostatic and magnetic fields external to the satellite 
until it is identified as either having hit some part of the 
satellite thus representing a returning fraction of the total 
current or else as having escaped, thus representing a current 
fraction to the environment. 


Particle pushing is done by solving the Lorentz force equa- 
tion using the leap-frog scheme of Boris. See Reference 4. 
Boris 1 procedure has been expanded to permit re-centering 
of the equations in time. Thus the timestep may be dynamic- 
ally increased or decreased as appropriate to keep parti- 
cles moving at velocities commensurate with the NASCAP grid 
in which they are located. 

k k , 

Tracking of particles which pass out of the highest NASCAP 
grid continues using a monopole electric field, the magni- 
tude of which is obtained from the total charge on the 
satellite. 
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9.4 


DISCRETE PARTICLE EMISSION ANGLES AND ENERGIES 


The present version of NASCAP offers the user two 
choices for emission angle selection. One choice is the uni- 
form distribution, a special case of which results in each 
particle representing the same solid angle fraction of the 
current. The other choice is a cosine 0 distribution in 
which a disproportionate number of particles are emitted at 
angles "close" to the axis of the hemispherical cap (Z-axis 
of the emitter coordinate system ) . Two choices for the energy 
spectrum of the beam are also provided. Either choice re- 
sults in an approximate representation of a mono-energetic 
peak in the energy spectrum — the difference between the two 
choices being the mathematical form of the approximation 
function. The emission angles and energy distribution func- 
tions available are listed below. (Any angular dependence 
may be combined with any energy dependence.) 

9.4.1 Uniform Angular Current Density Dependence 


For each of n discrete energies, n, n Q particles are 

£ (p y 

emitted. The initial emission velocity direction vector of 
each particle - (measured in the emitter coordinate system ) is 

V.. = (sin0 . cosd . ) i + (sin0 . sintf) . ) i + (cos0.)k (9.1) 

l j J i J i j 

where 


cj). = — (i - 1/2) 
r i n, 

9 


l X / 2 j w * • / n • 


e 1 = ^ (j - i/2) 

J u Q 


j = 1, 2, 


n r 


For the special choice of n^ = n^ each particle repre- 
sents the same solid angle fraction of the emitted current. 
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9.4.2 Cosine 0 Angular Current Density Dependence 

For each of n discrete energies, n, n Q particles are 

£ <p o 

emitted. The initial velocity direction vector of each parti- 
cle (measured in the emitter coordinate system ) is 


V . . = {sin0 . coscfi . ) i + (sin0 . sincj) . ) j + (.cos0 . ) k 
13 31 3 3 


where 


*i 


9 j 


~ 1 / 2 ) 


. -1 / (j-1/2) si-nA0 
sin [ — i 

0 


i = 1, 2 , . . . , n 


: = 


1 , 2 , 


(9.2) 


. . , n. 


9.4.3 Gaussian Energy Current Density Dependence 


A Gaussian function may be used as an approximation to 
a mono-energetic spectrum. The current density function for 
the Gaussian approximation is 


J(e) 



/2tt a 


exp 



where 1^ is the total beam emission current. 
0 


It is easy to show that 



ds 



Since a real current distribution is not defined for 
negative energies the Gaussian function is only an approxi- 
mate representation. However, one can show that 
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Lim J Ce) 
ff^-0 


I B 5 ~ £ o^ ' 


and 



de - I B 


for 


e >0 and — << 

O £ 

o 


1 . 


Thus the mono-energetic energy peak can be represented 
to any degree of accuracy desired simply by choosing °/^ Q small 
enough. It is worth noting that ^68 percent of the current 
falls in the range e = £ q + a and ^92 percent of the current 
falls in the range e = e q + 2a. 

NASCAP chooses the discrete energy representation of the 
Gaussian energy distribution as follows: 


e. = P -1 (j. + (1/2 - i) 

i \ n 


i = 1 , 2 , . . . , n 


where 


F(X) 


O 

-L 


/2tt o 


exp < - T 


1 I o 


2 V cr 


de 


(.9.3) 


This choice results in the following equality being 
satisfied : 


/ 


£ i+l 


J (e) de = — 
n 

£ 


(.9.4) 


where we define e =0 and e = ». 

o n £ +l 
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Thus each discrete energy e^ represents- the same frac- 
tion (l/n £ ) of the total emitted current. Furthermore, half 
of this fraction is a result of energy in the range e^_-^ < 
e £ e. and half is a result of energy in the range e^ <_ e < 


9.4.4 Rational Energy Current Density Dependence 


A rational function may be used as an approximation to 
a mono-energetic spectrum. The current density function used 
by NASCAP for the rational approximation is 


J(e) = 


y i 


B 


ir/2 + tan ^ y 


1 - 2 

e e 
o 


(■ 2 --S 

2 

1 + 

( e o sa ) 


(9.5) 


where I = total emitter beam current 
B 


y 


= { 7 ~- 


a /a and a = a/e 


This density function has the property that 


Lim J(e) -*■ I_ 5 (e - e •) 
a+0 B ° 


and 


J(e) de — I_ for all a/e < /2 
B ' o 

NASCAP chooses the discrete energy representation of the 
rational energy distribution according to Eq. (9.3) with the 
integrand for F(X) replaced by the rational density function 
divided by I_. This choice of energies also results in the 
satisfaction of the equality given by Eq. (9.4). 
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9 . 5 


NASCAP PARTICLE EMITTER ACTIVATION 


In order to activate low density particle emitters a 
special keyword file must be prepared. Each of up to five low 
density particle emitters is defined by the appearance of an 
emitter keyword definition sequence on this file. Each emitter 
keyword definition sequence consists of an "EMITER" keyword 
card which is optionally followed by one or more of the key- 
word cards described below. An emitter keyword definition 
sequence is terminated either by an "END" keyword card or else 
by another emitter keyword definition sequence. If more than 
one emitter is to be activated the associated keyword definition 
sequences must be consecutive. The last keyword sequence to 
appear on the emitter definition file must always be terminated 
by an "END" card. 

Once the special emitter keyword file has been prepared 
all that is necessary to enable the emitter (s) during the 
charging portion of subsequent TRILIN calls is to place a 
card on the NASCAP keyword input file {done using RDOPT from 
the runstream) which has the keyword "EMITTE" in columns 1-6 
and the FORTRAN unit number of the emitter keyword file in 
columns 29-30. If the emitter (s) have been activated by a 
previous RDOPT/TRILIN sequence then they may be deactivated by 
the keyword "NOEMIT" (appearing in the NASCAP keyword file) for 
subsequent TRILIN charging calculations. 
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9.6 NASCAP EMITTER SCRATCH FILES 

♦ 

The NASCAP emitter routines require the use of two 
scratch files . The NASCAP file IAUN is used as a temporary 
file for storing emitter fluxes to surface cells and other 
associated data. In general, the user need not concern him- 
self with this. The other file used is the NASCAP file IPART. 

It is the user's responsibility to make sure this file has 
been assigned for scratch usage. Since the file is used for 
particle trajectory plotting it is recommended that additional 
space beyond the computing system default for normal files 
be allocated. The default value of IPART is -28. (The nega- 
tive value is required to indicate that an old NASCAP routine, 
PARPLT , is not to be activated. The absolute value of IPART 
is used by the emitter routines, however.) On the Univac 1100/81 
machine it is recommended that the runstream include the 
card 

@ASG , T 28 , F///1000 
prior to the @XQT card. 
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9.7 EMITTER KEYWORD FILE PREPARATION 

At the beginning of each emitter keyword sequence 
(.appearance of the "EMITER" card in the emitter keyword file) 
all parameters which describe the properties of the emitter 
being activated assume a set of default values. The NASCAP 
user then has the option of changing any of these values to 
suit his requirements. This is accomplished by including one 
or more of the cards to be described below. The contents of 
each card consists of a mnemonic keyword, left- justified in 
card columns 1-6 and possibly the value of a data variable 
associated with the keyword. The type of data required (if 
any) is determined by applying the standard FORTRAN convention 
for variable types to the keyword. Values of variables must 
appear in columns 21—30 of the keyword card. The format used 
is 110 for INTEGER and F10.0 for REAL data. 
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9.8 


DESCRIPTION OF EMITTER PARAMETER OPTIONS BY KEYWORD 


9.8.1 

ICELL 

NPHIS 

NTHETS 

NENGS 

NSTEPS 

DTHETA 

JTYPE 


Emitter Dependent Parameters 


Index number of the surface cell at which an emitter 
is to be activated. The acceptable range is 1 < 

ICELL £ 1250. The default value is ICELL = 1. 

Number of discrete azimuthal angles at which parti- 
cles are emitted (Na in Eqs. (9.1) and (9.2)). Ac- 
ceptable range is NPHIS £ 1. Default is NPHIS = 1. 

Number of discrete polar angles at which particles 
are emitted (Ng in Eqs. (9.1) and (9,-2)). Acceptable 
range is NTHETS £ 1. Default is NTHETS = 1. 

Number of discrete energies at which particles are 
emitted (N £ in Eq. (9.3)). Acceptable range is 
1 £ NENGS £ 50. Default is NENGS = 1. 

Maximum allowable number of discrete time steps each 
of the NPHIS *NTHETS*NENGS emitted particles may be 
tracked before concluding a particle has escaped to 
the environment. Acceptable range is 1 £ NSTEPS £ 
2500. Default is NSTEPS = 500. (Only the first 
800 points of a particle's trajectory can appear in 
particle plots.) 

Polar angular width (in degrees) of the hemispherical 
cap used to model the emitter gun aperture (A0 in Eqs 
(9.1) and (9.2)). Acceptable range is 0.0° £ DTHETA 
£ 90.0°. Default is DTHETA - 1.0°. 

Index number of desired emitter gun current density 
function. The acceptable range is 1 £ JTYPE £ 4. 
Default is JTYPE = 1. The types are as follows: 

JTYPE 

1 Uniform angular, Gaussian energy distribution 

2 Cosine 9 angular, Gaussian energy distribution 

3 Uniform angular, Rational energy distribution 

4 Cosine 8 angular, Rational energy distribution 

Remark: For a detailed description of the current 

distribution see the section entitled "Discrete 
Particle Emission Angles and Energies". The central 
peak energy s 0 of any distribution is set by the key- 
word ENERGY and the width of the distribution in 
energy is specified by the keyword SIGMA. The 
Gaussian energy distribution has the attractive 
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property that a.68 percent of the current emitted is 
due to particles in the range of ENERGY + SIGMA. 

The NENGS discrete energies are selected so that an 
equal amount of the emitted current is represented 
by each energy. 

ENERGY Energy at which the central peak of the emitter cur- 
rent density function occurs. Acceptable range is 
ENERGY > 0.0. Default is ENERGY - 1000.0 eV. 


SIGMA 


ISPEC 


PHI 


Dispersion of the emitter current density function in 
eV (a in Eqs . (9.3) and (9.5)) . Acceptable range is 

SIGMA >0.0. Default is SIGMA = 0.1 eV. 

Emitter particle species type. The acceptable values 
are ISPEC = 1 for an electron emitter and ISPEC = 2 
for a proton emitter. Default is proton emitter 
(ISPEC = 2) . 

Azimuthal angle <j> of the emitter (in degrees) as 
measured in the spherical coordinate system of the 
cell at which the emitter is located. Acceptable 
range is 0.0° <_ PHI <_ +360.0°. Default is PHI = 
0 . 0 °.* 


THETA 


BEAMI 

VEDOWN 


Polar angle 0 of the emitter (in degrees) as measured 
in the spherical coordinate system of_ the cell at 
which the emitter is located. Acceptable range is 
-90.0° < THETA <_ +90.0°. Default is THETA = 0.0°.* 

Total emitter beam current in amps. Acceptable 
range is BEAMI > 0.0. Default is BEAMI = 10“6 amps. 

Tolerance limit for the maximum distance in code units 
which a particle is permitted to move at each timestep 
of particle tracking. If the limit is exceeded the 
timestep will be halved successively until the toler- 
ance limit is met. Acceptable range is VEDOWN > 0.0. 
Default value is VEDOWN = 0.3. (The limit is auto- 
matically scaled by the factor 2 iG- ^ where iG is the 
index of the ’grid in which the particle is located.) 


Note: Letting PHI and THETA assume their default values and 

setting DTHETA = 0.0 will produce an emitter which points 
along the normal to the surface cell ICELL upon which it is 
located . 
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VEUP 


(.VEDOWN /VEUP) is the tolerance limit for the minimum 
distance in code units which a particle must move at 
each timestep of particle tracking. If the limit is 
exceeded the timestep will be successively tripled 
until the limit is met. Default value is VEUP = 5.0. 
(If the default for VEUP and VEDOWN is used then the 
minimum distance a particle in grid 1 may move per 
timestep is 0.3/0. 5 - 0.06 inner grid units. As with 
VEDOWN, VEUP is also scaled by 2 iG-l.) 


Normally the boundary of a particle trajectory plot is 
adjusted to correspond to the boundary of the highest grid into 
which any of the NPHIS*NTHETS*NENGS emitted particles fell during 
tracking. Unfortunately particles circling in small magnetic 
fields tend to make very large loops, i.e., into grid' 12 or higher. 
Thus the automatic plot boundary selection frequently obscures de- 
tails of individual particles returning to the object in grid 1. 
(The object silhouette is only plotted if the boundary of the 
plot corresponds to grid 6 or lower and fine details are usually 
only apparent if the boundary of the plot is set to grid 4 or 
lower.) In addition there may be some cases in which it is 
desirable to have the same set of particle trajectories plot- 
ted more than once, each time using a different maximum grid 
boundary for the plot. The following three keywords provide 
the user with the capability to deal with these problems and 
special requirements. (Note that these three keywords only 
have an effect if JCYCEM >_ 1 — see Section 9.8.2.) Also note 
that these three keywords correspond to REAL data variables. 


ALPHA Grid number to which the boundary of the particle 

trajectory plots corresponds to. (Note that all 
grids which are higher than 2 but fall inside the 
plot edge are automatically outlined for reference 
on the trajectory plots.) Acceptable range is 
integral ALPHA >0.0, i.e., ALPHA = 0.0, 1.0, 2.0, 
... . If ALPHA = 0.0 then the boundary of the 

plot will be automatically adjusted to correspond 
to the highest grid into which any of the NPHIS* 
NTHETS*NENGS particles fell during tracking. De- 
fault value is ALPHA = 0.0. 
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BETA 


Number of different plots of same set of trajectories 
desired (only has an effect if ALPHA >_ 1.0) . Each 
trajectory plot is made using a different grid bound- 
ary (see GAMMA) . The acceptable range is integral 
BETA such that 1.0 <_ BETA <_ 4 . 0 . The default is 
BETA =1.0. 

GAMMA Increment factor for grid boundaries of successive 

trajectory plots. If BETA > 1.0 then the same tra- 
jectories will be drawn on BETA different plots. 

The^ grid boundaries are incremented for each succes- 
sive plot as follows : 

IG = ALPHA + (1-1) * GAMMA for 1=1, ..., BETA 

Acceptable range is integral GAMMA >_ 0 . 0 . Default 
is GAMMA =0.0. 

Note that the highest grid value selected for plotting 
boundaries by ALPHA, BETA, and GAMMA should not exceed LIMGRD 
if useful plots are to be obtained. 
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* 

9*8.2 Emitter Independent Parameters 


SCALEV Scale factor for the tolerance limit VEDOWN for parti- 
cle tracking in grids higher than 2. (Note that 
VEDOWN/VEUP is also affected by this.) Default is 
SCALEV = 1.0. Acceptable range is SCALEV > 0.0. 

LIMGRD Highest grid in which particle tracking is permitted. 

If particle passes outside this grid it will be as- 
sumed to have escaped to infinity. Default value is 
LIMGRD = 6. Acceptable range is LIMGRD 1. (Note 
that particles which exit from the highest NASCAP 
computational grid** are tracked using a monopole 
potential . ) 

PRFLUX If this keyword appears in the emitter file then a 

listing of all non-zero flux contributions to surface 
cells due to each particle emitter will be printed at 
the completion of each particle tracking cycle. (It 
is recommended that this card be included.) 

IPRNT If this card is included than an optional one-line 

summary is printed for each particle at the comple- 
tion of tracking. The information appearing in this 
line is as follows: 

IPHI index number of discrete azimuthal angle at 
which the particle was emitted. 

ITH index number of discrete polar at which the 

particle was emitted. 

IEK index number of discrete energy at which the 

particle was emitted. 

VINIT initial code velocity with which particle was 
emitted (in inner grid units/timestep) . 

VIN initial velocity with which particle was 

emitted (in meters/second) . 


If any of these keyword cards appear in any one of the one 
or more emitter keyword definition sequences on the NASCAP 
keyword input file then the parameter value will apply to 
all activated emitters. 

■k k 

The current version of NASCAP restricts particle tracking 
using explicitly calculated potentials from TRILIN to the first 
two grids even if potentials were computed in more than this 
number of grids . 
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J CLAST 


Index number of volume cell which particle 
was in at step just before hitting the satel- 
lite. If the trajectory was incomplete this 
will be 0. (Note the volume index is not 
necessarily the same as the surface cell in- 
dex! ) 

PXYZ Potential (in volts) at the particle posi- 
tion at the last timestep completed prior to 
hitting the satellite or abandoning tracking. 

IR Index number of grid in which the particle 

was in at the last timestep completed prior 
to hitting the satellite or abandoning 
tracking . 

ISTP Number of discrete -steps which this particle 
was tracked for before it hit the satellite 
or tracking was abandoned. 

Note in cases where the message *EMISSI0N SUPRESSED* 
appears prior to the trajectory end summary the 
values for JCLAST, PXYZ, IR, and ISTP are not cor- 
rect. Also if the warning that 50 emission suppres- 
sions have taken place is printed then values of 
these variables for trajectories which follow may 
not be correct. Setting the value of the INTEGER 
field on this card to 1 or 2 will produce a myriad 
of additional diagnostic output and is not recom- 
mended unless only one or two particles are being 
tracked. In general more than sufficient informa- 
tion can be obtained by leaving the data field 
blank (0 ) . 

CYMULT This parameter sets the "cyclotron" time limit. If 
a particle passes out of the highest grid* then this 
parameter essentially specifies the number of revolu- 
tions in the magnetic field which particles are per- 
mitted to make before concluding that they have 
escaped. That is to say tracking beyond the time 
at which a particle exits from the highest grid is 
done for at most (CYMULT* (2irm/ (eB) ) additional seconds 
before concluding that the particle has escaped. In 
special cases where there is no magnetic field pres- 
ent this parameter has a different use. Let T be the 
time which a particle requires to pass out of the 
highest grid*. If it passes out of this grid and 


* 

loc. cit. 
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there is no magnetic field then it will be tracked 
for a maximum of T*CYMULT additional seconds before 
concluding it has escaped. Acceptable range is 
0.0 < CYMOLT £ 10.0. The default is CYMULT = 1.0. 
(In cases of no S-field it will probably be neces- 
sary to set CYMULT >5.0 if premature tracking term- 
ination is to be avoided.) 

JCYCEM 'Emitter particle trajectory plot flag. JCYCEM = 0 
implies no plots. JCYCEM ■> 1 implies that particle 
trajectory plots will be produced for all emitters 
every JCYCEM time cycles performed by TRILIN, begin- 
ning with cycle 1. Default is JCYCEM = 0. 

IPLTYP Particle trajectory plot type selection. IPLTYP = 0 
implies separate sets of trajectory projection views 
for each emitter. IPLTYP = 1 implies all emitter 
trajectory plots combined in a single set of three 
two-dimensional projection views . IPLTYP only has 
an effect if JCYCEM > 1. Default is IPLTYP = 1. 

COMMEN This card is for programmer convenience in identify- 
ing the purpose of keyword cards on the keyword in- 
put file. If this card appears, columns 7-72 may 
contain any comment the user desires. It is echoed 
back on the keyword listing at execution time but 
is otherwise ignored . 


9.9 KEYWORD CARD PROCESSING ERRORS 

Following each emitter keyword definition sequence NASCAP 
performs an extensive set of checks to assure that an acceptable 
set of parameters has been defined. If the set is found un- 
acceptable then diagnostic information is printed for each 
offending parameter and program execution is aborted. If no 
errors are detected than a summary of all parameters for each 
activated emitter is printed on the NASCAP output listing. 


9.10 SPECIAL RESTRICTIONS 

1. If more than one emitter is simultaneously activated no 

two emitters may have the same underlying conductor. 

No error check for this condition is explicitly made 
so that results may be unpredictable if this constraint 
is violated. 
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9.11 

1 . 


2 . 


ADDITIONAL OPERATION NOTES 

Most of the printed output produced by the emitter 
routines is self-explanatory. The particle emitter 
energies printed are the actual energies with which, 
each of NPHIS *NTHETS particles are emitted at. The 
particle emitter energy weights printed on the output 
currently have no function within the program and 
should be ignored by the user. 

In general, the quantities printed following the 
heading ***** SUMMARY OF PARTICLE TRACKING FOR EMITTER 
XX - AT - CYCLE - XXX **** have been derived from a 
consideration of all NPHIS*NTHETS*NENGS particles which 
were emitted. 

A number of diagnostic **WAKNING** statements may be 
printed during particle tracking. These warnings are 
only intended to be informative about how the tracking 
algorithm is progressing and should only be cause for 
concern if something "bad" happens afterward such as an 
***ERROR*** message and/or execution abort call to 
RETRNO before execution of the emitter routines is 
completed during the charging cycle of TRILIN. 
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10. AMBIENT CHARGE DENSITY 


A local net charge density can develop in a neutral 
ambient plasma near an object because of the unequal attrac- 
tion for the electrons and ions in the plasma. Our treatment 
assumes that the net charge density developed in the ambient 
plasma is linear in the local potential. Thus Poisson's equa- 
tion, 


e o V <j> = -p 


( 10 . 1 ) 


becomes 


where 1 is an effective Debye length. Laframboise and Parker 
have shown ^ 5 ^ that for a Maxwellian plasma, the correct linear 
term is generated using 


I* ♦' 


(. 10 . 2 ) 



2 - 1/2 

where = (4irn e /kT ) ' is the true Debye length, n is 
the plasma density far from the object, and T^ and T^ are the 
electron and ion temperatures, respectively. We have used 
these approximations to develop an ambient space charge cor- 
rection for NASCAP . 

We must find a variational functional which is appro- 
priate for the governing equation of the system, Eq. (10.2) . 
The surface terms enter the equations as previously derived, 
and we must reconsider only the volume term. The correct free 
volume functional is 
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c- !> 



f V = | e o (V<j>).(V<f>) + 


(10.4) 


The requirement of v = 0 yields Eq. (10.2). In a finite element 
formulation we have 


f V - f {§ « 0 * i * jwl * m3 + I ( 72 ) n1n3 


Extremizing with respect to <j> yields 


dV 


t e i j + X ± j = 0 


(10.5) 


( 10 . 6 ) 


where 


£ 


ij 



VN 1 • VN D dV 


(10.7) 



Matrix elements L. . .have been calculated for each of the 
5 distinct standard cell types - ; these new elements are given at 
the end of Chapter 10. The matrix elements e. . are of course 
identical to the corresponding terms derived in the absence* of 
space charge. The iterative solution of the potential equa- 
tions proceeds exactly as before, but with the elements 
[e^j + A^j] entering wherever only occurred previously. 

We have also recast the boundary conditions for the 
problem to be compatible with the presence of the ambient 
sheath. When IOUTER = 2 is selected, the routine SETALL now 
sets the outer boundary to 
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where R q is an average radius for the object, and Q is the 
charge on the object. Equation (10.9) represents the correct 
analytical solution of Eq. (10.2) for a spherical object. 

The above ideas have been implemented by the following 
changes or additions to the indicated NASCAP routines: 

SUBROUTINE FLXDEF: Additions to calculate the ef- 

fective Debye length, A. 

SUBROUTINE SETALL: Additions to calculate R , the 

average object radius, and to set 
the boundary potential according 
to Eq. (10.9) . 

SUBROUTINE SCLIN : New routine to calculate the matrix 

sum [e^j + A^j] , given A. 

The routine SCLIN is called once per timestep by POTENT, and 
the resulting matrix sum is stored in COMMON/WGTS/ to be 
accessed by routines ECUBE and SQCWGT during the potential 
iterations . 

Minor changes in routines SQCWGT, ECUBE, POTENT, and 
all the input routines were required as well. A new keyword, 
"ISC", has been introduced; the above procedures are invoked 
by inputting a value of 1 for ISC. 
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10.1 RESULTS 


Sample results are shown ,in Figures 10.1 - 10.4. Fig- 
ures 10.1 - 10.3 show potential contours around an aluminum cube 
of edge 2 m as the Debye length decreases from 33' m to 1 m. The 
cube is fixed at -100 V while the boundary is grounded in these 
calculations, and the potential contours run from -5 to -100 
volts in all figures. The development of the ambient sheath is 
clearly reflected in the changing potential contours. Figure 
10.4 shows the change of capacitance of a 1 m aluminum cube in 
unbounded space as the Debye length changes. The ability of 
the object to store charge is dramatically increased by the 
presence of the ambient sheath, as expected. 

Notice that the amount of computation required for each 
potential iteration is essentially unchanged when the linear 
ambient correction is included, since the matrix sum 
is calculated before the time-consuming iteration procedure is 
entered. Since the matrices A^j are more positive definite than 
the corresponding e. ., the convergence is rapid, and fewer itera- 
tions are required. The overall computing requirement is there- 
fore usually reduced when the ambient terms are included, so 
that we have been able to extend the range of applicability of 
NASCAP without a concomitant increase in running time- 
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POTENTIAL CONTOURS ALONG THE X-Y PLANE OF 2 • = 



Figure 10.1. Potential contours around a 2 m aluminum cube, 
including the effect of ambient charge density. 
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POTENTIAL CONTOURS ALONG THE X-Y PLANE OF Z = P 

ZnlN = -.10000"- 03 ZflAX - -.120*5+00 A 2 = .50000+01 



Figure 10.2. Potential contours around a 2 m aluminum cube, 
including the effect of ambient charge density. 
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POTENTIAL CONTOURS ALONG THE X-Y PLANE OF Z = 9 

31IN = 10000*03 ZMAX * -.S7J34-04 iZ = .50000+01 



Figure 10.3. Potential contours around a 2 m aluminum cube, 
including the effect of ambient charge density. 
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farads) 



Figure 10.4. Capacitance of a 1 m aluminum cube versus Debye 
length. 
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MATRIX ELEMENTS FOR SQUARE OF POTENTIAL- 

(Format) 

Description 
Standard Orientation 

Potential Function = N 1 ^ 

i 

Matrix, L^ ^ : J'&tt d 2 = ^ L ij^i^j 

ij 

Point Index 

1 
2 

3 

4 

5 

6 

7 

8 


Cube Corner 

0 0 0 
10 0 
0 10 
110 
0 0 1 
10 1 
Oil 
111 
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Standard Cell 0 
Empty trilinear cube 
Orientation: Arbitrary 


Potential Function: 


i 

N 1 


1 

(1- 

>1 

1 

r~ I 

2 

(1- 

•z) .(1-y): 

3 

(1- 

■x) y (1-z 

4 

(1- 

■z) yx 

5 

z(l-y) (1-x 

6 

x(l-y) z 

7 

zy (1-x) 

8 

xyz 

(0) 



ij 



1/27 



1/54 

1/27 


1/54 

1/108 

1/27 

1/108 

i/54 

1/54 

1/54 

1/108 

1/108 

1/108 

1/54 

1/216 

1/108 

1/216 

1/54 

1/216 

1/108 

1/108 





1/27 




1/216 

1/27 



1/108 

1/54 

1/27 


1/108 

1/54 

1/108 

1/27 

1/54 

1/108 

1/54 

1/54 


1/27 


195 




Standard Cell 1 

Half-Empty Wedge 

1 < x + y < 2 
0 < z < 1 



Orientation: Right angle along 

line 7-8 


Potential Function: 


1 
1 

2 

3 

4 

5 

6 

7 

8 


0 

(1-y) (1-z) 
(1-x) (1-z) 
(x+y-1) (1-z) 
0 

(1-y) z 
(1-z) z 
(x+y-1) z 



0 


0 

1/3 6 





0 

1/72 

1/3 6 




0 

1/72 

1/72 

1/9 



0 

0 

0 

0 

0 


0 

1/72 . 

1/144 

1/144 

0 

1/3 6 

0 

1/144 

1/144 

1/144 

0 

1/72 

0 

1/144 

1/144 

1/18 

0 

1/7 2 


1/3 6 

1/72 1/9 
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Standard Cell 2 


Cube with diagonal line on one face 
Orientation: Line from 2 to 3 



Potential Function: 
i N 1 

1 (1-x-y) (l-z) (1-x-y) 

2 [x0 (1-x-y) + (1-y) 0 (x+y-1) ] (l-z) 

3 [y0 (1-x-y) + (1-x) 0 (x+y-1) ] (l-z) 

4 (x+y-1) (l-z) 0 (x+y-1) 

5 (1-x) ( 1-y) z 

6 x (1-y) z 

7 ‘ (l-x)yz 

8 xyz 


As a simplification, these cells have been treated as 

if they were type 0 cells. Previous experience with matrix 

2 

elements of j V<J> | indicated that the calculated potentials 
were insensitive to such an approximation. 

as defined previously. 
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Standard Cell 3 
Tetrahedron 

2<x + y+ z<3 



Orientation: Empty corner at 

point 8 

Potential Function: 


1 
1 

2 

3 

4 

5 

6 

7 

8 


N 
0 
0 
0 

1-z 

0 

i-y 

1-x 

x+y+z-2 


l.' 3 > 

13 


0 

0 0 


0 

0 

0 





0 

0 

0 

1/60 




0 

0 

0 

0 

0 



0 

0 

0 

1/120 

0 

1/60 


0 

0 

0 

1/120 

0 

1/120 

1/60 

0 

0 

0 

1/120 

0 

1/120 

1/120 


1/60 
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Standard Cell 4 


Truncated Cube 

Orientation: 000 corner (point 1) 

missing 



As a simplification, these cells are treated as if they 
were type 0 cells. A value for the potential at point 1 is 
generated by extrapolation from the potentials at points 2 
through 8 , the matrix elements for type 0 cells is applied 
(multiplied by a factor of 5 / 6 to account for the reduced 
volume) , and the results are interpolated back to the seven 
real points. The results of this procedure are equivalent to 
introducing a matrix for type 4 cells given by 


L (4) =* (5/6) T t L (0) T 


where T is a matrix which performs the extrapolation, and L 

<V rs» 

is the matrix defined previously for type 0 cells. 


The extrapolation scheme used is 


= P + (P - P ) 
1 m m a 


where 

P m = -j (P 2 + P 3 + P^) = potential at midpoint of 

triangular face 

P a ~ T ^ P 4 + P 6 + P 7 + = avera 9 e potential behind 

the triangular face 

This particular scheme was chosen since it yields elements a 

( 4 ) 

matrxx L . . which is positive definite, and whose elements are 
1 J . (4) 

all non-negative. The exact must satisfy these require- 

ments . 
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Thus the matrix T is given by 


T 

11 

= 0 

T ii = lr 

2 < i < 8 



T 

12 

= T 13 = 

T 15 = 2/3 




'S' 

i — 1 

EH 

=. t = 
16 

T 17 + T 18 

= - 1/4 



All 

other T 

H- 

l_i. 

1! 

o 

. 




The 

resulting is 

• 

• 



(4) 
L . 
13 






. 00000 






.00000 

.06516 





.00000 

.04201 

.06516 




.00000 

.01157 

.01157 

.02894 



.00000 

.04201 

.04201 

.00000 .06516 



.00000 

.01157 

.00000 

.00579 .01157 

.02894 


. 00000 

.00000 

.01157 

.00579 .01157 

.00579 

.02894 

.00000 

.00129 

• .00129 

.01447 .00129 

.01447 

.01447 


.03086 
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11 . SURFACE CONDUCTIVITY 


The first part of this chapter (.Section 11.1) is- a paper 
by Rotenberg, Mandell, and Parks reproduced verbatim. This 
gives an analytical discussion of bulk and surface conductivity. 

The second section (11.2) shows how the analytical model 
is incorporated into the NASCAP code. 
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11.1 EFFECTS OF BULK AND SURFACE CONDUCTIVITY OF THE POTENTIAL 
DEVELOPED BY DIELECTRICS EXPOSED TO ELECTRON BEAMS* 


** 


Manuel Rotenberg , Myron J. Mandell and Donald E. 

Systems, Science and Software 
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ABSTRACT 


The charging and discharging of a dielectric material 
which has bulk and surface conductivities is discussed. Two 
model problems are solved. In the first problem, a semi- 
infinite dielectric plane, attached to an infinite grounded 
conducting substrate and exposed to a monoenergetic electron 
beam, is analyzed. Bulk and surface conductivities, and 
secondary emission characteristics are taken into account as 
parameters . In the second problem the dielectric is charged 
but the electron beam is shut off so only the bulk and sur- 
face conductivities enter the calculation. The principal re- 
sult of the latter calculation is to show that steep tangen- 
tial gradients develop in the presence of surface conductivity 

during decay, and that for asymptotic times the temporal be- 

- 1/2 

havior, for a fixed position, is proportional to t ' , 

rather than exponential as expected in the presence of bulk 
conductivity . 
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1 . INTRODUCTION 


This study of charging and discharging of dielectric 
materials is motivated by the importance that has been attached 
in recent years to the effects produced by charging of satel- 
lite surfaces exposed to the magnetospheric environment. Under 
substorm conditions surfaces are exposed to electrons with tem- 
peratures of the order of 10 keV and thermal plasma currents 
of the order of 1 nanoampere per square centimeter. Thus, on 
a time scale of a few minutes, large potential differences 
(^10 keV) may develop between the surface of' the thin dielectric 
(^10 cm) and an underlying conductor, producing electric 
fields of the order of 10 6 volts/cm. At such field levels, the 
occurrence of electrical discharges, with concomitant RF noise 
and material degradation effects, is expected. Such expecta- 
tions are confirmed by both laboratory and space data. One 
suggested means for mitigating these effects is to develop 

dielectric materials with sufficient (bulk or surface) con- 
ductivity to suppress the development of excessive electric 
fields. 

Laboratory experiments investigating charging and dis- 
charging phenomena typically involve exposure in high vacuum 
of thin dielectric films mounted on conducting substrates to 

several kilovolt electron beams carrying currents in the 
2 

nanoamp/cm range. The dielectric materials and configura- 
tions tested include kapton, teflon, nonconductive paints, 
solar cell arrays, second surface mirrors and optical solar 
reflectors. The electric fields which develop in the bulk of 
the dielectric and near the interface between dielectric films 
and metal substrates are a determining factor in the occur- 
rence of discharges. The fields in turn are determined by 
the electrical characteristics of the material, in particular 
bulk and surface conductivities and secondary emission yields. 
However, it is very difficult to measure by conventional means 
the very low conductivities typical of spacecraft materials. 
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In the present study, we consider the role of con- 
ductivity and secondary emission in establishing the temporal 
and spatial dependence of potential difference V{x,t) between 
the dielectric surface and the .underlying conductor in condi- 
tions which simulate the space environment- In particular the 
analysis describes the spatial and temporal variation of 
potential near the edge which separates the dielectric coated 
and exposed metallic surfaces, and the manner in which it de- 
pends on material properties . 

When the charging current is turned off, the analysis 
of the discharging behavior provides a basis for measurement 
of bulk and surface conductivities. In addition, it will be 
shown that the resistance of the dielectric-conductor edge 
interface, if it is significant, can be measured indirectly. 
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2 . DIELECTRIC CHARGING 


Consider a thin dielectric on a grounded conducting 
plane (Figure 1) . The dielectric thickness is small compared 
with any lateral dimension. Thus the electric field perpen- 
dicular to the plane is V/d. Parallel to the surface of the 
dielectric the field is 3V/3x. (We suppose no y-variation. ) 

Let j be the net current density (incident minus secondary 
nst 

emission) impinging upon the dielectric in a perpendicular 
direction. 

The incident current is assumed to be uniform over the 
surface of the dielectric. The continuity equation reads 


3u _ . 1 V + 1 3 2 V 

3t -*net p d w 2 


( 1 ) 


where 

a is the charge density 
p is the bulk resistivity 

w is the surface resistivity 

d is the thickness of the dielectric 

If the dielectric is semi-infinite in extent with its edge at 
x = 0, and is initially uncharged, the boundary conditions for 
(1) are 

V(x, 0) = 0; V ( 0 , t) = 0; 3V(«,t)/3x = 0 . (2) 

In typical cases, both the intensity of incident current 
and its energy (thus secondary yield) are functions of the sur- 
face potential. Secondary yield is a non-linear function of in- 

g 

cident energy, therefore to cast Eq. (1) into a tractable form, 
we assume 3 net ( x ) varies linearly near its zero value (see Fig- 
ure 2) as a function of local potential 

j net (x) = f[V 0 " V(x)] (3) 
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where V q is the equilibrium potential for a perfect dielectric 
( w = p = co) of infinite extent, and use Gauss' Law in the form 

a = V/d (4) 


where « is the dielectric constant. Using Eqs . (3) and (4), 

(1) becomes 


3V 

at 


-kV + D 




(5) 


where k is defined as the inverse time constant 


k = 


d 

K£ o 




( 6 ) 


and 


D 


__d 

K S CO 

o 


(7) 


is defined as the diffusion constant. It is convenient to 
write (5) in the dimensionless form 


aw 

8x 


-W + 



+ 1 


W(£,0) = 0; W(0,x) = 0; 3W(«>,t)/a£ = 0 


(8a) 

(8b) 


where 


x = kt 

S = (k/D) 1/2 x 

W = v/(rv o ) 

F = Dyu/k = ypd/Cypd + 1) 


( 9 ) 



Equation (8) is solved in a straightforward manner by Laplace 
transforms : 


W(£,t) 



erf 



+ 


1 

2 



erfc 





er fc 




( 10 ) 


This function is plotted in Figure 3. 

The asymptotic solutions to Eq. (8) are 


W(C,°°) -► 1 “ e 


and 


(ID 


W(°°,t) -> 1 - e T (12) 

so it might be thought that a useful approximate solution to 
(8) is 


W(?,x) - (1 - e“ ? ) (1 - e T ‘) (13) 

Comparison with the exact solution shows that Eq. (13) is a 
poor approximation. In particular, (13) fails to exhibit the 
rapid rise to a constant value as a function of E, for x < 1. 
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3. DIELECTRIC DISCHARGING 


Assume the incident current is shut off after the di- 
electric reaches its equilibrium potential distribution. This 
time is defined as t — 0 for the discharging problem. The 
governing equation is (5) with y set to zero. For convenience, 
the normalized variables used in the charging problem will be 
retained: 


9W 

9x 


-s W + 


o 



The boundary conditions are now 


(14a) 


W(£,0) = 1 - e a? ; W-(0 ,t) 


0; 9W(co,t)/9£ = 0 


(14b) 


where s^ = k /k, k = l/(pKe ). The distance scale a has 

been introduced into the initial condition to compensate for 

beam deflection as the dielectric is being charged. It will 

2 

be seen that for x >> £ , the results are independent of a. 
The solution to Eq. (14) is 


W(5,t) = e 


-s x 
o 


erf 


2x 


1/2 


i 2 - 

1 ax 

" 2 e 


+a$ 


e a ^erfc ( ax 


1/2 


2x 


1/2 


- e erfc ax 


1/2 


2x 


1/2 


(15) 


1/2 

which is plotted in Figure 4. If £ << x ' it is seen that 
the W is a linear function of the appropriate limit of (15) 
yields 

1/2 


-s x 

W(£,x) = e ° 5/(-irx) 


(16a) 
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or, in terms of unsealed, variables. 


V(x,t) = V exp [-t/(pK£ )] (C w) 1//2 (16b) 

° ° (TTt) X/ 

where C = Ke Q /d is the capacity per unit area of the di- 
electric . 

Equations (15) and (16) , together with Figure 4, con- 
stitute our main result. The exponential factor is expected 
for the discharge through a dielectric of non-zero bulk con- 
ductivity when the surfaces are uniformly charged initially. 
When there is surface conductivity, however, the charge be- 
comes non-uniform toward the edge of the dielectric and re- 
sults, for a fixed x, in the inverse square root behavior in 
time. This non-uniform charge distribution is of practical 
importance since it produces large potential gradients tan- 
gential to the dielectric surface, the existence of which 
could lead to sporadic surface discharges. 
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4 . SURFACE CURRENT 


If the conducting substrate of the dielectric is 
grounded through a current measuring device, the surface 
current at x = 0 is measured. (It will be shown subsequently 
that the bulk current is small compared to the surface cur- 
rent for a reasonable choice of parameters . ) The surface 
current is given by 




1 9V 
a) 9x 



exp 


-t/(pKS Q ) 



(17) 


5 

The time constant pK£ Q is typically mlO seconds so the ex- 
ponential factor can be replaced by unity. 


211 



5 . BULK CURRENT 


The bulk current density is given by 


= h V(5 - T) 

If V is taken to be its maximum value, rV Q , then 

I /I < area 1 / tt(o \ 1//2 t l/2 

x bulk / surf perimeter p \Ke^6.) 


area , -2 ,1/2 , , s 

t — x 10 t (meters) 

perimeter sec 


(19) 


where it has been assumed that 


p = 10 


16 


SI -meters 


,-4 

w = 10" Si/a 


d = 10 meters 
13 


If the dielectric sample is taken to be a square 10 cm on a 

-3 1/2 

side, the ratio of bulk-to-surface current is about 10 t ' . 

4 sec 

Thus observations could be made for about 10 seconds before 
the bulk current became an appreciable fraction of the sur- 
face current. The area-to-perimeter ratio could be improved 
by using a smaller sample, but if this were carried too far, 
the one-dimensional analysis breaks down. 
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6 . DISCUSSION 


The time and space dependence of two problems have been 
analyzed; the charging of a semi-infinite dielectric film by 
a uniform beam, and the subsequent discharge of the film after 
the beam has been turned off. 

The basic assumptions in the charging problem are that 

t 

the secondary emission is a linear function of the beam energy, 
and that the electrical properties of the dielectric can be 
clearly divided into bulk and surface sensitivities. There is 
a further assumption that the beam current is independent of 
time and spatial coordinates. 

The discharge problem retains only the assumptions re- 
garding the resistivities and for this reason may be more re- 
liable. Two independent measurements can be made to test the 
assumptions: potential measurements near the dielectric edge 

to test Eq. (16b) , and the current measurements from the 
conducting plate to ground to test Eq. (17) . 

Potential measurements near the edge of the dielectric 
may not extrapolate to zero at'-x = 0 as demanded by Eq. (16b) ; 
this would be an indication that the edge of the dielectric 
presents a non-negligible resistance to the surface current. 
Thus suppose that the voltage at x = 0 is found by extrapola- 
tion from interior measurements to be V^(t), and that the 
total current leaking to the conducting plate across the edge 
of the dielectric is I (t) . Then the resistance of the edge 

is R = V"(t)/I(t). If the edge of the dielectric is A meters 
e o 

long then the relevant intrinsic quantity is the edge conduc- 
tivity per unit length 

g = (1/A) (I/V^) mho/meter (20) 

The independence of g over time is a necessary behavior (but 
not sufficient) for the edge, bulk and surface resistances 
to be independent of potential. 
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Jsurf 
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x=0 


Figure 1. Definitions of the coordinate system and the various 
currents. The net incident current is j ne .(- = jj_ nc 
- jsec* The dielectric has a thickness d and a di- 
electric constant tc. 



Figure 2. A typical secondary current characteristic. V 0 is 

the point at which j net = 0. The slope of the curve 
at this point is defined as y. 
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Plot of Eq. (10) , the solution to th 
line showing the value of the normal 
value at £ = 00 was found numerically 
from the edge of the dielectric; the 
dielectric as a fraction of the ulti 


lelectnc c 
d potential 
The absciss 
dinate is t 
e voltage. 


roblem. Tne 
percent of its 
e scaled dista 
acre across the 









v/rv, 



a£ 

Figure 4. The solution to Eq. (14); the normalized potential as a function of 

normalized distance from the edge of the dielectric. The parameter is 
normalized time. The abscissa is the scaled distance from the edge of 
the dielectric; the ordinate is the voltage across the dielectric as a 
fraction of the initial voltage. 


11.2 SURFACE CONDUCTIVITY 

NASCAP treats surface conductivity as a redistribution 
of charge among the nodes of a surface cell due to bilinearly 
interpolated surface fields in that cell. The resultant 
matrix elements are shown on the following page. For an in- 
sulating material, surface conductivity is activated by 
specifying a positive surface resistivity value (fi/e) as 
material property 14 . This also produces conductivity at a 
dielectric-metal edge. A negative value indicates surface 
conductivity is to be ignored. The ignoring feature is im- 
portant, since there is a limit of 9537 matrix elements in the 
circuit model. 


SELECTED SURFACE 

RESISTIVITIES (fl/n) 

Teflon (TFE) 

10 13 

Teflon (FEP) 

>10 16 

Kapton 

10 16 

Mylar 

10 16 

RTVxxx 

10 12 

Epoxy 

1-t 

H 

o 

rH 

Glass (quartz) 

io 19 


218 



Surface Conductivity m.e.'s 


Units are 
a 

s 

£ o 

A 

Square 


Rectangle 


Right 

Triangle 


Eq. 

Triangle 


/e A , where 
s o 

surface conductivity (mhos) 
— 12 

8.85 x 10 farads/meter 


mesh size (XMESH) 



'll 


a 


12 


13 


-2/3 

^14 ~ ^ / ® 
1/3 


/2 



= -7 /2/12 
= / 2/12 
= /2/4 
= /2/4 


= a 13 " 1/2 

" a 33 = ~ 1/2 
= 0 


0 22 a 33 4/3/3 

" a 13 = a 23 = 2/3/1 
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12. DISCHARGE, PATCH, SPACE CHARGE DENSITY PLOTS 


12.1 THREE FURTHER MODIFICATIONS 

These three code modifications proved to be simpler 
tasks than many of the previously described extensions. Yet 
they significantly expand NASCAP capabilities. 

The DISCHARGE facility identifies sites of spacecraft 
discharge. It involves a new keyword (read by RDOPT) and two 
new material properties (read by OBJDEF) . 

PATCHR and PATCHW are two new object types to be used 
in defining complicated objects. They are very much like the 
RECTANGLE and PLEDGE objects. The difference is the way they 
interact with the hidden line HIDCEL routine. 

Space charge density plots can be obtained using the 
keyword SHEATH (read by RDOPT) . 

12.2 NASCAP DISCHARGE ANALYSIS 

A discharge facility has now been implemented in NASCAP. 
To invoke it, specify the RDOPT option DISCHArge, with the re- 
laxation factor DCHGF in columns 21-30. The breakdown poten- 
tials are specified as material properties 15 and 16 (both 
positive) : 

Property 15: Maximum absolute potential attainable by 

material . 

Property 16 : Maximum potential of a dielectric surface 

relative to its underlying conductor, or 
a non-grounded conducting surface relative 
to spacecraft ground. 

Analysis is performed immediately preceding return to 
TRILIN from LIMCEL. In the case of a plasma discharge (prop- 
erty 15 exceeded) charge is removed from the lowest numbered 
non-fixed conductor, as well as charge transfer from dielectric 
surface to underlying conductor. If property 16 is exceeded. 
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the obvious charge transfer takes place. The amount of charge 
transfer is governed by the "relaxation factor" DCHGF, with 
DCHGF = 1. indicating total short , and DCHGF = 0. indicating 
minimum charge transfer. 

Example 1 . A grounded aluminum cube, partially covered with 
teflon, is exposed to a charging environment. The teflon sur- 
face is set to discharge to its conducting substrate with a 
50 percent relaxation factor when differentially charged to 
1 kV. The final portion of LIMCEL printout for times tep 4 
(126 seconds to 286 seconds) is shown in Figure 12.1. During 
this time the teflon surface charged from v-89 0 volts to 
^-1800 volts, exceeding the discharge threshold. Discharges 
were registered on all the teflon surface cells, and charge 
was transported from the surface to the conductor. At the 
beginning of timestep 5 the teflon surface was seen to have a 
potential of ^-500 volts. 

Example. 2 . The partially covered aluminum cube of example 1 
was allowed to float. The aluminum surfaces were set to dis- 
charge to test tank ground at a potential of 5 kV, with a 
70 percent 'relaxation factor. During timestep 13 the aluminum 
charged from -4700 volts to -5900 volts. (The teflon surfaces 
remained 500 volts negative relative to the aluminum.) As 
shown in Figure 12.2, a discharge was registered at the lowest 
numbered aluminum surface cell. At the beginning of timestep 
14 the aluminum cube was at 'v-lSOO volts. 
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PREDICTED NODE POTENTIALS 


NODE 

V-PRED 

30 

1 

-1 b64. 36 

-469210. ID 

2 

-1732.91 

-928232.57 

7 

-J 

-1664.38 

-469208 .42 

-~r 

-1732.92 

-923231 .45 

5 

-1767.27 

-922546.40 

6 

-1732.9 1 

-92S202 .55 

7 

-1 S06 . 6 £ 

-1 535623.16 

s 

-1732.92 

-928203.65 

9 

-1806,69 

-1 S3 5 6 2 1 . 30 

10 

-1804 .34 

-1333437.08 

1 1 

-1664.38 

-469203. 42 

12 

-1732.92 

-9232-31 .47 

13 

-1767.27 

-922546.43 

14 

■-1732.92 

-928230.65 

15 

-1806.69 

-1 835621 .30 

16 

-1SC4 .84 

-1833437.05 

17 

-1767.23 

-922544 .52 

18 

-1804 . 3 4 

-1833435.39 

1 9 

-18C3.35 

-1374053.16 


V-OLD p T L I S T 

-512.38 001010131001 
-349. SO 001011131001 
-512.89 001012111001 
-849.60 001012121001 
-363.15 001012131 OQL 
-349.60 001110131001 
-889.09 001111130001 
-849 .61 00 1 1 121 1 1001 
-989.09 001112120001 
-883.50 001112130001 
-312.39 001210111001 
-849.60 001210121001 
-363.15 001210131001 
-549.61 001211111001 
-389.09 001211120001 
-885.50 001211130001 
-363.16 031212111001 
-888.50 001212120001 
-S37.96 031212130331 


SUM ( D Q ) - -2.21 256995 + 07 


S U'M ( C S ^ C V ) - -1 . 4 37 92 35 6 +04 

DOC - 2.211 13212+07 

V-C - 0.G0C00QC0 

DISCHARGE ANALYSIS 50. UX RELAXATION 

CELL 4 AT -i .73+03 VOLTS DISCHARGE TO UNDERLYING CONOUCTOR 

CELL S AT -1.73+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 10 AT -1.78+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 11 AT -1.78+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 15 AT -1.73+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 16 AT -1.76+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 18 AT -1.73+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 20 AT -1.7S+C3 VOLTS DISCHARGE TC UNDERLYING CONDUCTOR 

CELL 21 AT -1.78+33 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 22 AT -1.60+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 23 AT -1 .80+ 03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 

CELL 24 AT -1.80+03 VOLTS DISCHARGE TO UNDERLYING CONDUCTOR 


Figure 12.1. Discharge printout - grounded aluminum cube. 


•^asr 


222 



PREDICTED NODE POTENTIALS 


NODE 

V-PRED 

DO 

V-CLD 

PTLIST 

1 

-6465.78 

-197.85 

-5250.02 

001010131001 

2 

-6449.86 

-397.47 

-5232.71 

001011131001 

3 

-6465.88 

-197.85 

-525C.10 

0010121 11Q01 

4 

-6449.89 

-397.48 

-5232.73 

001012121001 

5 

-6423.79 

-400.08 

-5206.93 

001012131001 

6 

-6449.86 

-397.47 

-5232.71 

001110131001 

7 

-6433.80 

-798 .36 

-5215.96 

0C1111 130001 

8 

-6449.91 

-397.49 

-5232.74 

001112111001 

9 

-6433.31 

-798 .40 

-5215.97 

001112120001 

10 

-6414.39 

-602 .40 

-5196 .76 

0 G 1 1 12130001 

11 

-6465 .88 

-197.86 

- 5 2 5 G . 1 0 

001210111001 

12 

-6449.89 

-397.49 

-5232.73 

001210121001 

13 

-6423.79 

-4 00 .-08 

-5206.93 

001210131001 

1 4 

-6449 .91 

-397.49 

-5232.74 

OC 1211111 OCi 

15 

-6433.81 

-798 .40 

-5215.97 

001211120001 

1 6 

-64 14 .39 

-602 .-40 

-5196.76 

0G1 21 1 130001 

17 

-6423.85 

-400.12 

-52C6.97 

001212111001 

IS 

-64 14 .40 

-802.45 

-5196.77 

001 212120001 

19 

-6401 .58 

-603.99 

-5184 .30 

0Q1 212130001 


SUM ( OQ.l -9 . 5651 I 5 60+0 3 
SUM { CS^DV ) - -i .96546336+04 
DOC = -1.60014203+04 
V-C = -5.91873303+03 

DISCHARGE ANALYSIS 70. OS RELAXATION 

CELL 1 AT -5.9Z+03 VOLTS DISCHARGE TO PLASMA OR TANK GROUND 

Figure 12.2. Discharge printout — floating aluminum cube. 
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12.3 PATCH SURFACES 

The object definition routines INPUT , RECTAN, and WEDGE 
have been rewritten to accept two new object types: "PATCHR" 

and "PATCHW" . These two objects are identical to the previous 
objects "RECTAN" and "WEDGE", respectively, except that no 
A2-type (shadowing) surfaces are defined for PATCHR and PATCHW. 
The new objects can be used to change the surface cell composi- 
tion of selected areas on large, previously defined surfaces 
when the outline of the entire object is not to be altered. 

Before the new keywords were implemented, the list of 
A2 surfaces was becoming unreasonably large for objects with 
complicated surface composition patterns, such as SCATHA. The 
changes described above allow the A2 surface list to remain 
sensibly related to the geometrical complexity of an object, 
rather than to the complexity of the surface composition. 

12.4 SPACE CHARGE DENSITY 

The NASCAP space charge density plotting subroutine 
RHOSHE may be invoked by the keyword SHEATH. It obtains to 
first order the space charge due to low energy photo- and 
secondary-electrons emitted from surface cells (not, as yet 
booms) . The contour levels are in units of code units of charge 
per cubic mesh unit. The effect of space charge on the space- 
craft potential (in volts) is comparable to the maximum contour 
level . 

3 

To convert a contour level, z, to MKS units (coul/m ) , 

2 

multiply by £ q /A , where A is the mesh size in meters: 

q £ Z q 

p (coul/m ) = — y— 'v 10 z 
A 
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3 +6 

To obtain electrons per cm , divide by a further x 10 

, , , 3, 10~ 6 £ o ' . 0 4 

p (el/cm ) = — r 2 ^ 1° 2 

g e 

Figures 12.3 and 12.4 show 'space charge contours for 
the small SCATHA model in an extreme environment. The first 
picture is for an uncharged satellite, and shows that the 
largest space charge is in the cavity, where it has a value 
of z = -.29 = -5x10 11 coul/m 3 = 300 el/cm 3 . The second is 
for the spacecraft at -675 volts, and shows no space charge 
except in the cavity, where it attains a peak value half that 
of the first figure. 
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13 . SHADOWING CALCULATIONS 


Systems, Science and Software performed a calculation 
of percentage solar illumination for SCATHA experiments- The 
resulting tables show experiment illumination for the various 
solar orientations that SCATHA will assume. Detailed results 
are given in Systems, Science and Software Report SSS-R-78-3658 , 
May 1978, "SCATHA Experiment Shadowing Study". 

The source code and relocatables and data files were 
delivered to NASA-LeRC. This section documents those items. 


13.1 CODES AND DATA FILES 

A Univac 1100/80 16-track program file tape has been 
delivered to NASA-LeRC. This tape contains nine program files, 
any of which could be reloaded to mass storage using the @COPY,G 
command. Of the nine files on the tape the ones which are sig- 
nificant to NASA-LeRC are probably the first five. It is sug- 
gested that NASA-LeRC retain copies of all nine files, however, 
for possible future reference. The contents of the first five 
files are as follows: 

1. SHADOW This file contains all of the source and re- 

locatable elements for the shadowing program 
as well as an absolute element generated on 
the Univac 1100/80 at S 3 . The absolute element 
contains S 3 plot package routines, however, and 
it is recommended that NASA-LeRC remap the pro- 
gram. (See next section.) 

2. STARTRUN This is a symbolic element file which contains 

all of the canned runstreams used to generate 
the 33 SCATHA shadow tables. 

3. OBJSCA This symbolic element file contains polygon 

definition keyword cards for each of the ob- 
jects on SCATHA which were significant for 
the shadow tables requested. 

4. AlOBJ This symbolic element file contains combina- 

tions of elements from OBJSCA used to make up 
the complete shadowed object polygon sets. 
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5. A20BJ This symbolic element file contains combina- 

tions of elements from OBJSCA used to make up 
the complete shadowing object polygon sets. 

3 

The S 1100/80 run which copied the nine files to tape 
is reproduced in Figure 13.1 and should assist NASA-LeRC in 
making their own copy of the tape. 


1 

4 

5 
A 


i o 
1 1 
is 

1 3 

14 

15 
IS 
1 7 
19 
i 9 

90 


9M?R PLEBBE MOUNT. -MEM UNLR.BELEII S'.RVE TRPE RB PEP POST 
9BCS*TJ TRPP. IJ9 
9PFUTNTI TRPE. 

9MRPk TRPE. 

3MRPU TRPE. 

9PP1.IT NO TRPP, 

9 COPY . RM OHROOW * TRPP 
9 COPY * CM BTRPTPUN. TRPP 
9 COPY* RM OB. Jo CH * ThPP 
9 COPY * RM R 1 OB J * TRPE 
9 COPY * RM RPOBJ * TRPP 
9 COPY * RM 3RMTEC * TRPE 
9 COPY * RM RVEPRRE * TRPE 
9C0PY * RM NEMMRT N. TRPE 
9 COPY * RM SCRo HBR01.I * TRPE 
3MRPf«- TRPP. 

9MRPK’ TRPE. 

9 PPM TNTi TRPE. 

9 PPM I Nil TRPP. 

9PPEE TRPE. 


Figure 13.1. Generation of shadow program tape sent to NASA- 
LeRC. S3 tape number was N1976 (unlabeled) . 
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13.2 REMAPPING AT NASA-LeRC 


3 

The shadow programs utilize the S routines 

S3DATE, 

S3ETIM, 

S3TICK, 

S3WARN. 

If equivalent routines are not available at NASA/LeRC then the 
following dummy subroutines could probably be implemented with 
the only loss in flexibility being the capability of writing 
restart dumps periodically or writing a restart dump and grace- 
fully terminating before exceeding the CPU time limit. 

SUBROUTINE S3 DATE (DATE, TIME) 

DATE = 0 
TIME = 0 
RETURN 
END 

SUBROUTINE S3ETIM ( ZTIME) 

ZTIME = 0 

RETURN 

END 

SUBROUTINE S3TICK (TIME) 

TIME = 0.0 

RETURN 

END 

SUBROUTINE S3 WARN (* , LIMTIM, LIMPAG) 

RETURN 

END 
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Once these dummy routines have been written and com- 
piled they should be mapped with the relocatables in the 
SHADOW file using the following procedure: 

Assume that the relocatable for the NASA/LeRC version of 
the four routines discussed above are already in TPF$ . Then 

§COPY , R SHADOW. 

@MAP 

@COPY , A NAME $, SHADOW. ABS 
@PACK SHADOW. 


If undefined references to the routines SETUPV, SETUAV , DXDYV 
GRID1V, TYPEV, and possibly DTLINE appear they should be ig- 
nored. 



13.3 DESCRIPTION OF PROGRAM RUNSTREAM INPUT 


Immediately following the @XQT SHADOW. ABS control card 

the program expects to find the following in the runstream: 

Card 1: 6 character file identifier selected by the user. 

Appears in card col. 1-6. 

Card 2,3: Two 78 character run description lines — any text 

chosen by the user. The text should be centered on 
these cards for neat appearance on the output list- 
ing. 

The remaining data cards in the runstream consist of 

one or more sets of the following card sequences: 

A. Task command keyword card. (Card col. 1-7 left- justified . ) 

B. SHDATA NAMELIST data cards to select options for the 
keyword command specified by A. 

The available keyword commands and their functions are: 

NEWFILE: A brand-new shadow table (or segment of) is to be 

generated. At completion of this command the table 
will be written onto FORTRAN logical file number 
LUNNEW. 

RESTART: Resume calculation using an old shadow table re- 

siding on FORTRAN logical file number LUMOLD which 
was not completely filled during a previous run 
due to excessive execution time. 

EXTEND: Extend an old shadowing table. The old table re- 

sides on FORTRAN logical unit LUNOLD and the hew 
table containing the old entries plus the new seg- 
ment calculated is written onto file LUNNEW. 

MERGE : 
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Shadow tables on files LUNOLD and LUNNEW are merged 
into a new file which is written onto LUNMRG. No 
actual shadow calculations result from this command. 



PRINTO : 


Produce a line printer listing of the shadow table 
residing on file LUNOLD . 

PRINTN: Produce a line printer listing of the shadow table 

residing in LUNNEW. 

PRINTM: Produce a line printer listing of the shadow table 

on LUNMRG. 

INSERT: Insert entries into a shadow table "by hand". (Fol- 

lowing the S HD AT A namelist another namelist FIXUP 
is also expected. This namelist contains the one- 
dimensional array SHTAB of shadow entries and NFILE — 
the number of entries in the table. The dimension 
of SHTAB is SHTAB (NPHI ,NTHETA) . If NPHI = 360, 

NTHETA = 26, for example, and we wish to set 
SHTAB (231,14) = 0.31 then the FIXUP namelist would 
be 

$FIXUP SHTAB ( 5271) =0 . 31 , $END , 
for example. This command assumes the table to be 
added to already exists on LUNOLD. The new table 
will be written to LUNNEW. LUNMRG is used as a 
scratch file. 

END: All tasks for the run are complete. Terminate execu- 

tion. 

The variables which may appear in the SHDATA namelist, 

their functions, and default values are as follows: 

IUNIT1: Logical unit number upon which the shadowed object 

A1 polygon definition keyword cards reside. De- 
fault is IUNIT1 = 10. 

IUNIT2: Logical unit number upon which the shadowing object 

A2 polygon definition keyword cards reside. De- 
fault is IUNIT2 = 11. 

LUNOLD: Logical unit number of old shadow table file. 

LUNMRG: Logical unit number of merged shadow table. 
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LUNNEW: 

ICMPR: 


APHI: 


BPHI : 


NPHI: 


ATHETA: 

BTHETA: 

NTHETA: 


JP 

NP 

JT 

NT 


Logical unit number of new shadow table file . 

Flag set as follows: 

0: shadow table in full 2-D format. 

1 : shadow table in compressed format with zero 

entries removed. 

Currently ICMPR = 1 has not been fully implemented. 
Therefore the user should not attempt to input a 
value for ICMPR. The default is, of course, 

ICMPR = 0. 

Lower limit for azimuthal angle <> in table (in 
degrees). Default is APHI = 1.0. 

Upper limit for azimuthal angle <j> in table (in 
degrees) . Default is BPHI = 360.0. 

Number of azimuthal angles (rows) in table. 

Note: This parameter should, in general, remain 

the same for keyword command/SHDATA namelist se- 
quences which appear in the runstream after the 
first and in subsequent runs which involve the same 
table. Default is NPHI - 360. (This is the maxi- 
mum allowable value.) 

•Lower limit for polar angle 0 in table (in degrees) 
Default is ATHETA = 80.0. 

Upper limit for polar angle 0 in table (in degrees) 
Default is BTHETA = 100.0. 

Number of polar angles (columns) in table. 

Note: Remarks about NPHI apply here also. De- 

fault is NTHETA = 21. Maximum allowable value is 
NTHETA = 26. 

Starting index for azimuthal angles.. 

Ending index for azimuthal angles . 

Starting index for polar angles . 

Ending index for polar angles. 
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If a new table segment is being generated then the range 
of the table filled is 

SHTAB { I , J) for I = JP, . .., NP 

and 

J = JT , . . . , NT 


The angle corresponding to I, J is 


<f> ± = DPHI* ( I— 1) + APHI 

0 . = DTHETA* ( J-l) + ATHETA 
1 

where 

DPHI = (BPHI - APHI) / (NPHI - 1) 

DTHETA = (BTHETA - ATHETA) / (NTHETA - 1) . 


' INONES : If INONES ^ 0 set all entries in table in range 

specified by JP, NP, JT, NT to 1.0. (No actual 
shadow calculation performed. Facilitates inser- 
tion by hand of large areas of the table which are 
known before-hand to be total shadow.) Default is 
INONES = 0. 


LIST: If LIST f 0 then a complete listing of the table on 

file LUNNEW is produced following completion of 
shadowing calculations . 
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13.4 OBJECT DEFINITION FILES 


If a shadowing calculation is to be performed then the 
Al and A2 polygon definition files must be available to the 
program. The data cards on these files must be in the follow- 
ing format: 

First Card: *A1DEF or *A2DEF 

( polygon keyword definition cards for each 
object in the shadowed (shadowing) polygon 
set. 

Last Card: *ENDA1 or *ENDA2 

For the SCATHA satellite a special file (OBJSCA) has been 
prepared which contains polygon definitions for each of the 
SCATHA objects deemed significant to the shadowing study. The 
name of each element in this file, in general, corresponds to 
the associated SCATHA object which is defined by the element. 

The individual objects for a particular shadowing table are 
selected from the file OBJSCA and combined (using the @'ED pro- 
cessor of the Univac 1100) into shadowing object ( A2 ) and 
shadowed object (Al) polygon sets. A number of these sets can 
be found in the files AlOBJ and A20BJ. In many cases, however, 
the elements needed from OBJSCA were simply combined using @ED 
within the runstream used to generate the shadow table. 
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13.5 RUNS TRE AMS 


A special symbolic element file called STARTRUN has been 
provided with the program. This file contains all of the canned 
runstreams used to generate the 33 SCATHA shadow tables. 

Many of the tables were generated segment-by-segment 
using several consecutive ©START commands using elements from 
the STARTRUN file. In most of these cases there is an ad- 
ditional ©COPY operation required between ©STARTs to make the 
output file from a previous run available to the next. 

As a simple example we shall consider the generation 
of the shadow table for the ML12-6 TQCM aperture. This shadow 
table was generated in four segments. All segments were 
generated within the same runstream, however, so only a single 
batch run was required. The runstream for this table is con- 
tained in STARTRUN . ML6AP , a listing of which is reproduced in 
Figure 13.2. The appropriate shadowing object sets for ML12-6 
are as follows : 


SHADOWING OBJECT SET AND APPROXIMATE SUN ANGLE RANGES: 


Shadowing Objects 


Approximate Sun Angles 


1. SC2-1 (B) . boom and sheath 
electric field sensor 

2. SCll-l(B) boom and magnetic 
field sensor 

3. Main satellite body: The 

ML12-6 radiator surface and 
aperture are not visible to 
the sun for the range 


4. ML12-6 door: Test Calculations 

showed that shadowing was im- 
portant only at sun angles for 
which 0 > 100°. Thus the door 
was not included in the 
shadowing object set for the 
ML12-6 surfaces. 


23° 

< 

<f>- < 

48° 

89° 

<r 

e < 

97° 

188° 

< 

* < 

201 

90° 


0 < 

95° 


0° < <j> < 22° 
202° < <*> < 360° 
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3 ED S.MLGRP 
RERD-ONLY HOSE 
ED S14. 00-04^20- 
EDIT 
OsLNP* 


.8t39-sOy> 


4 

5 

8 
“7 : 
I 

8 

g 

10 

It 

12 

13 

14 

15 
18 
17 
13 
13 
20 
21 
22 
i— _/ 

24 

25 
28 
27 
23 
23 
3 0 
81 
io 

JU 

30 

34 

35 
33 




■5RUN PGS y 1 1 073" 0 0 y STEEN—P * • 

3ES+U . E IG-NRMEy N 

SCRTHR SHflDOU 

TRELE FDP ML12 

CDNTRfllNRTIEJN 

EXPERIMENT =:XML6RP> 

3USE R1 >-RlDEJ 
3IJSE R2>R20EJ 
5USE □ * QEJ8CR 
3USE 2 * SCRSHRDGU 
3USE N* NEUMRIN- 
SRSG* T 10, 

SRSGyT 11. 

3RSG>PUP XE3LB. 

3RSG> PUP XNEU. 

3U8E 20* XOLD 
SUSE 22>XNEU 
3ED □. BLRNKy 10. 

I ♦RIDEF 
RDD Q.MLia-GIN 
I +ENDR1 
LNP* 

5ED D.ELRNKyll. 

I +R2DEF 
RDD O.SC2-1BOGM 
RDD 0. SC2— 1SPHEPE 
I *ENDR2 
LMP l 

5R8GvT 12. 

5EB Q. BLANK v 12. 

I ♦R2DEF 
RDD □. SCI 1—1 
I ♦ENDRE 
LNP? 

®RSG*PUP XML6RP. 

-5XQT N. 

XML6RP 


;600 




ML 12-6 RPERTUPE SHRDDUED B' 


;C2— 1=.E> AND SCll-l :E> 


Figure 13.2 


Runstream used to generate shadowing table for 
ML12-6 aperture. 
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4 Li 
41 : 
4c! 
40: 
44: 

45 

46 

47 

48 

49 

5 0 
51 
58 

50 

54 

55 


58: 

59 

60 
61 
68 
60 

64 

65 

66 
t y i* 


NEUFILE 

S3HDRTR 

JUT I NE=i 00, LIMTI M-3 0 , 

JP=23,NP=49» 

JT =8 , NT =8 0 > 

send 

EXTEND 

SSHDRTR i_UN0LD==82 , LUNNEU=2 0* IUMI T2= 
JP=184,NP=801, 

JT=1 , NT=2i, 

SEND 

EXTEND 

S3HBRTR LUN0LD=2 0 , LUNNEU=£2 , 
IMOMES^l , 

JP=8 08 v NP=36 0 , JT= 1 , NT=2 1 , 

SEND 

EXTEND 

S3HBRTR LUNDLD=22,LUNNEU=£0, 
INDHES=1, 

JP=1 , NP=£2> JT=1 v NT=81 » 

SEND 

END 

SPNDvGED 
5C0PY 20, XML6RP 
9FP.EE Xf1L6fiP 
3X0T RVEPRGE. 

SfiVG PNOR!t=0. 0,-0.3746, 0.9878, SEND 
®ES*U.BIG~NRNE,N 
CHRBOU 


6 y 

70 

7 1 


TRELE 
END 
BOX o: 
9FIN 


3CRNJ72 


EOF :72 
0 : EX I 

NO CORRECTIONS RPPLIED. 



The polygon definition for the ML12-6 aperture is found in 
OBJSCA. ML12-6IN. The required Al polygon definition file 
is prepared in lines 18-22 of the runstream using the editor. 
The polygon definition for the SC2-l(Bj boom and sphere are 
found in OBJSCA. SC2-1B00M and OBJSCA. SC2-1SPHERE, respectively. 
Since these two objects shadow over a different range than the 
SC11-1 they are placed on their own file at lines 23-28. The 
definition for SC11-1 is found in OBJSCA. SC11-1 and is placed 
on another file by lines 29-34. Line 35 assigns a permanent 
file onto which the final shadow table is to be copied. At 
line 37 program execution begins. (The absolute program ele- 
ment was assumed to be in file NEWMAIN — at NASA/LeRC this 
will be in the file SHADOW instead.) After reading the 
identifiers in lines 37-39 the first segment of the table in- 
volving shadowing by SC2-1(B) is generated by lines 40-45. 

The second segment involving shadowing by SC11-1 is 
generated by lines 46-50. Note that the output file from lines 
40-45 becomes the input file for lines 46-50. Finally the 
main satellite body shadowing effect is inserted in two seg- 
ments by lines 51-55 and 56-60. (Note the exchange of input 
and output shadow table files again.) The END command word 
at line 61 terminates execution and the completed shadow table 
is copied onto the. save file at line 63. Line 65-66 uses the 
completed table as input to compute and print out the inte- 
gral averages requested by Dave Hall for this experiment. If 
a printed copy of the- final shadow table is desired it could 
be obtained by setting LIST = 1 at line 57 or by making a sub- 
sequent PRINTO run with the file XML6AP assigned to unit 20 — 
the LUNOLD file. 
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14. MATCHG - A MATERIAL CHARGING CODE 


14 . 1 INTRODUCTION 

MATCHG is a zero dimensional code that charges a surface 
using the same material formulations found in the NASCAP code. 
The surface potential, V, is just 

V = Q/C , 

where .Q is the charge per unit area and C is the capacitance 
per unit area 



The code models either a monotonic electron beam or a spherical 
probe approximation to a Maxwellian flux of ions and electrons. 
The initial version contains no conductivity effects. 

Output includes tables and plots of electron backscatter 
and secondary emission versus energy and the surface potential 
versus time. 

Section 14.2 describes how to use the code and Section 
14.3 contains two sample runs. 
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14.2 USE OF THE CODE 


Material data is read off unit 8. The format is that 
used in the NASCAP object definition file, i.e., a material 
name followed by three lines of properties. The code allows 
properties to be changed interactively so that the material 
properties file need not be modified. MATCHG is designed for 
an 80-line terminal. On EXEC 8, this can be set by @@TTY W,80. 
The code will prompt the user. Below are listed appropriate 
responses to questions. In general YES is the only affirmative 
response recognized and a carriage return (cr) will suffice 
for a negative response, but NO., or anything but YES, will also 
be interpreted as a negative. 

To terminate execution reply to the MATERIAL prompt with 

STOP . 

Prompts and acceptable responses are listed below. 


PROMPT 


RESPONSES 


EMISSION FORMULATION 


MATERIAL 


PRINT 


{ cr, leaves it unchanged 
ANGLE-, regular formulation 
NORMAL , normal emission formula 
only 


cr, leaves it unchanged, if first 
call stops the program 
STOP, stops the code 


Provided in 
MATCHG . DATA 


ALUMINUM j 
TEFLON ( 

KAPTON ( 

SI02 ) 

MGO ' 

YES, lists material properties 
( must be called if properties 
are to be changed) 


NO 

cr 


No list 
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PROMPT 


RESPONSES 


CHANGE ANY MATERIAL 
PROPERTIES 


[property] 


FLUX TYPE 


NEW PARAMETERS? 


[flux parameter] 


GENERATE A TABLE? 


PLOT ( ) ? 


CHARGE? 


FULL PRINT? 


cr 

NO 

YES 


! cr — leave value unchanged 
[value] — new value 

( cr — remain unchanged (not first 
time) 

1 — monoenergetic electron beam 

2 — Maxwellian with electrons 

and protons 

f cr 

) NO 
( YES 

( cr — leave value unchanged 
) [value] — new value 


makes table of electron 
secondary, electron back- 
scatter, and proton generated 
electron secondary yields 
versus incident energy for 
normally incident particles 


makes plot 


/ cr 
NO 

YES - 


cr 

NO 

YES - 


cr 

NO 

YES — charges sample 

cr — prints first and last cycle 
NO current balances 

YES — prints every cycle potentials 
and every fifth cycle cur- 
rent balances 
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PROMPT 

CHARGING POSITIVE — 
STOP? 


PRINT A CHARGING TABLE? 


PLOT? 


244 



14.3 SAMPLE RUNS 


The first sample run is for ALUMINUM and shows most of 
the possible types of output. The second run is for KAPTON 
and shows how little output can be produced and still produce 
meaningful parameter variations. 
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93TTY 1 . 1 , 133 
-^COMPLETE 
> 3ASG., ft MPRP 
READY 

>?*USE 3,MPRP 
PEADY 

'•5COPY, A MCHG.ABS 

FURPUR 27R3-1 E35 SL73P.1 07x13/73 09:21:05 

1 APS 
>.?XQT 


WELCOME TO MATCHG - A MATERIAL CHARGING PROGRAM 


EMISSION FOPMULATION IS ANGLE 
SECONDARY EMISSION FOPMULftT ION'» NORMAL 
FMI^ION FORMULATION IS NDRMAL 
MATERIAL'. ALUMINUM 
PRINT >7ES 


original 

0f POOR 


page m 
QUALITY 


PREPROCESSING OF MATERIAL PROPERTIES 
MATERIAL, l: ALUMIN 



PPOPERTY 

INPUT VALUE 

CODE VALUE 

1 

DIELECTRIC CONSTANT 

l . nn+fio 

'NONE} 

1 . 00+00 

(NONE) 

2 

THICKNESS 

1. 00-03 

METERS 

1.00-03 

MESH 

3 

CONDUCTIVITY 

- 1 . 00+00 

MHO/M 

- 1 . 00+00 

MHO/M 

4 

ATOMIC NUMBER 

1. 30+01 

<NQNE} 

1.30+01 

■.NONE} 

5 

DELTA MAX '■COEFF 

9.70-01 

<N0NE3 

7.21+00 

* NONE} 

6 

E-MAX >DEPTH*+— 1 

3. 00-01 

KEV 

1.73-02 

ANG-01 

» 

PAN 6 E 

2.60+02 

ANG. 

3.33+02 

ANG. 

3 

EXPONENT > PANGE 

1 .30+00 

(.NONE} 

4. 15+02 

ANG. 

Q 

RANGE > EXPONENT 

2.40+02 

ANG. 

1.30+00 

(NONE) 

ib 

EXPONENT 

1.73+00 

(NONE) 

1.73+00 

(NONE} 

li 

YIELD FOP 1KEU PROTONS 

1.36+00 

'■NONE} 

1 . 36+00 

(NONE'* 

12 

MAX DE/DX FOP PROTONS 

4. 0 0+ 01 

KEV 

4. 00+01 

KEV 

13 

PHOTOCURPENT 

4. nn-05 

A/M ♦♦2 

4. 00-05 

A/M++2 

14 SURFACE RESISTIVITY 

CHANGE ANY MATERIAL PROPERTIES) YES 
DIELECTRIC CONSTANT >2.0 

0 . 00 

OHMS 

0. 00 

V-?/Q 


THICKNESS > 1 . E— 2 

CONDUCTIVITY >2 

ATOMIC NUMBER > 

DELTA MAX >CDEFF > 

E-MAX 7DEPTH++-1 > 

RANGE > 

EXPONENT > RANGE > 

RANGE > EXPONENT 
EXPONENT 7 

YIELD FOP 1 KEV PPOTONS > 

MAX DE/DX FOP PROTONS > 

PHOTOCURPENT 

SUPFACE RESISTIVITY > 


MATERIAL i: ALUMIN 


1 

PROPEPTY 

DIELECTRIC CONSTANT 

INPUT vrluE 
2 , 00+00 (NONE} 

CODE VALUE 
2. 00+00 (NONEV 

2 

THICKNESS 

1 . 00-02 

METERS 

1 . 00-02 

MESH 

3 

CONDUCTIVITY 

2 . 00+00 

MHO'M 

2 . 00+00 

MhO/M 

4 

ATOMIC NUMBEP 

1 . 30+01 

(NONE} 

1.30+01 

(NONE) 

5 

DELTA MAX ;-COEFF 

9.70-01 

iNONE) 

7.21+00 

CNONE'' 

6 

E-MAX 7 DEPTH+* — 1 

3. 00 -ni 

KEV 

1.73-02 

ftNG-01 

7 

RANGE 

2.60+02 

ANG. 

3. 33+02 

ANG. 

B 

EXPONENT 7 PANGE 

1 . 30+00 

CNONE) 

4. 15+02 

ANG. 

Cl 

RANGE 7 EXPONENT 

2.40+02 

ANG. 

1. 30+00 

(NONE) 

ib 

EXPONENT 

1.73+00 

(NONE) 

1.73+00 

•NONE'* 

n 

YIELD FOP 1KEV PROTONS 

1 . 36+00 

■.NONE'* 

1.36+00 

(NONE) 

12 

MAX BE/DX FOR PPOTONS 

4. 00+01 

KEV 

4. 00+01 

KEV 

13 

. PHOTOCURPENT 

4. 00-05 

A-T1++2 

4. 00-05 

8' M ++2 

14 

SURFACE RESISTIVITY 

0. no 

OHMS 

0. 00 

V— S/'*? 

ENTER FLUX TYPE 1=BEAM 2=MAX«ELLIAN >1 

BEAM VOLTAGE = 2. 0+01 KEV BEAM CURRENT = 1 . 

NEW PARAMETERS 7 >YE3 

0-05 AMPS 

/ M++2 



DERM VOLTAGE CKEV*>12 
BEAM CURRENT (AMP/M++2) 7 

BERM ■•'OLTRfiE = 1.2+01 KEV BEAM CURRENT = 1.0-05 AMPS 7 M++2 

GENERATE A TABLE ?>YES 



-RGYcKE'-O 

EL. SEC. 

EL. BkSCAT. 

PP-. SEC 

.100 

.300 

. 054- 

.429 

.2 00 

.333 

. 108 

.605 

. 300 

.370 

. 138 

.739 

.400 

.331 

. 159 

.352 

.50 0 

.355 

.175 

. 950 

• GOO 

-773 

. 187 

1.033 

.70 0 

.714 

. 197 

1.11.3 

.S00 

.353 

. 206 

1.193 

l . ono 

.577 

.219 

1.327 

2. 000 

.334 

.204 

1.832 

3. oon 

. 302 ' 

. 192 

2. 191 

4. 000 

.254 

. 182 

2.473 

5. 000 

.221 

. 174 

2.703 

6. 000 

.133 

. 167 

2. 397 

7. ono 

.130 

. 162 

3. 068 

3. 000 

• 136 

. 157 

3.206 

3. ono 

. .154 

. 154 

3. 331 

10. 000 

.144 

. 151 

3.441 

20. 000 

. 032 

. 139 

4. 055 

30. 000 

. 071 

. 137 

4.257 

40. 000 

. 058 

. 137 

4. 301 

50,000 

. 050 

. 137 

4.274 


PLOT SECDNBRRY ELECTRONS'? >YES 


'1 . 00 + 00 + 

I * 

I * 

I* 

3.5(1-01+ + 

r 

i ♦ 

i 

7.00- 03+ ♦ 

I + 

I + 

r+ + 

5.50- 01+ 

I 

I 

I 

4.00- 01+ + 

I 

I 

I * 

2.50- 01+ ♦ 

I + 

I 

I 

0.00 2.50+00 5. no+on 

YIELD •FRACTION’) Vi. ENERGY <"KEV'> 
PLOT BRCKSCATTERED ELECTRONS -?>Y ES 



1 

7.50+00 


1 . 0 0+0 1 


3. 00- 01+ 

T 

i 

I 

2.50-01+ 

I 

I + 

I + 

2 . 00 - 01 + + 


r + ♦ ♦ 

i ♦ 
i + 

1.50-01+ 

I ♦ 

I 

1 + 

1.00-01+ 

I 

I 

I 

5. 00-02+*- 
I 
I 
I 


0.00 2.50+00 5.00+00 

YIELD <FRRCTION'> '■>* . ENERGY CKEV.'' 
PLOT PROTON SECONDARIES ? >YES 


50+00 


1 . 00+01 
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6 . 00 + 00 + 

I 

I 

5. 00+00+ 

I 

I 

I 

4. 00+00+ 

I 

I 

I 

3. 00+00+ 

I 

I 

I 

a. 00+00+ 

i + 

i 

i+- 

l . oo+oo++ 



i 

'n.nn l.oo+oi a.on+m 3.00+01 4.00^01 5.00+01 

YIELD (FRACTION} VS. ENEPSY (KEV1 
CHARGE 8 >YES 
FULL PRINT ?>7ES 

CYCLE 1 TIME 0.00 SECONDS POTENTIAL 0.00 VOLTS 


INCIDENT ELECTRON CURRENT -1.00-05 

SECONDARY ELECTRONS 1.23-06 

BACKSCATTEPED ELECTRONS 1.46-06 
INCIDENT PPOTDN CURRENT 0.00 

SECONDARY ELECTRONS 0.00 


NET CURPENT -7.26-06 AMPSYM++2 


CYCLE 


TIME 

2.33—01 

SECONDS 

POTENTIAL 

-1.18+03 

VOLTS 

CYCLE 

3 

TIME 

5.86-01 

SECONDS 

POTENTIAL 

-2.34+03 

VOLTS 

CYCLE 

4 

TIME 

8.78-01 

SECONDS 

POTENTIAL 

-.3.47+03 

YOLTS 

CYCLE 

5 

TIME 

1.17+00 

SECONDS 

POTENTIAL 

-4.58+03 

'■'DLTS 

CYCLE 

6 

TIME 

1.46+00 

SECONDS 

POTENTIAL 

-5.64+03 

VOLTS 

INCIDENT 

ELECTRON 

CUPPENT 

- 1 . 

0 0-05 




SECONDARY ELECTRONS 1.31-06 

BACKSCATTEPED ELECTRONS 1.65-06 
INCIDENT PROTON CURRENT 0.00 

SECONDARY ELECTRONS 0. 00 


NET CURRENT -6.44-06 AMPS'M+*2 


CYCLE 

7 

TIME 

1 . 76+0 0 

SECONDS 

POTENTIAL 

—6. 66+03 

VOLTS 

CYCLE 

3 

TIME 

2.05+00 

SECONDS 

potential 

-7.62+03 

VOLTS 

CYCLE 

Cl 

TIME 

2.34+00 

SECONDS 

POTENTIAL 

-8. 51+03 

YOLTS 

CYCLE 

lb 

TIME 

2.63+00 

SECONDS 

POTENTIAL 

—9. 31+03 

VOLTS 

CYCLE 

1 1 

TIME 

2.93+00 

SECONDS 

POTENTIAL 

—9. 99+03 

YOLTS 

INCIDENT 

ELECTRON 

CUPPENT 

-1.00-05 




SECONDARY ELECTRONS 3.82-06 

BACKSCATTEPED ELECTRONS 2.04-06 
INCIDENT PPOTON CURRENT ■ 0.00 

SECONDARY ELECTRONS 0.00 


NET CUPPENT -4.14-06 AMPS/'M*+2 


CYCLE 

12 

TIME 

3.22+00 

SECONDS 

POTENTIAL 

-1 . 05+04 

CYCLE 

18 

TIME 

3.51+00 

SECONDS 

POTENTIAL 

-i . no-+ci4. 

CYCLE 

14" 

TIME 

3.S1+00 

SECONDS 

POTENTIAL 

-t . 11+04 

CYCLE 

15 

TIME 

4. 10+00 

SECONDS 

POTENTIAL 

-1.13+04 

CYCLE 

16 

TIME 

4.39+00 

SECONDS 

POTENTIAL 

-1 . 14+04 


INCIDENT ELECTRON CUPPENT -1.00-05 

SECONDARY ELECTRONS 7.57-06 

BACKSCATTEPED ELECTRONS 1.31-06 
INCIDENT PPOTON CURRENT 0.00 

SECONDARY ELECTRONS 0. 00 


VOLTS 

YOLTS 

YOLTS 

VOLTS 

VOLTS 


NET CURRENT -5.18-07 AMPS-'M++£ 

^oj-nri ■s'c-^rrwri-c* DnTCNTtfll 


CYCLE 

17 

TIME 

4.68+00 

SECONDS 

POTENTIAL 

-1.14+04 

VOLTS 

CYCLE 

13 

TIME 

4. 98+00 

SECONDS 

POTENTIAL 

-1 . 14+04 

YOLTS 

CYCLE 

19 

TIME 

5.27+00 

SECONDS 

POTENTIAL 

-1. 14+04 

YOLTS 

CYCLE 

20 

TIME 

5.56+00 

SECONDS 

POTENTIAL 

-1. 14+04 

VOLTS 

CYCLE 

21 

TIME 

5.36+00 

SECONDS 

POTENTIAL 

-1. 15+04 

YOLTS 


INCIDENT ELECTRON CUPPENT -1.00-05 

SECONDARY ELECTRONS 8.17-06 

BACKSCATTERED ELECTRONS 1.31-06 
INCIDENT PPOTON CURRENT 0.00 

SECONDARY ELECTRONS 0.00 


6 . 00+01 
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NET CUPRENT -1.52-nB AMPS^M++2 


CYCLE 

22 

TIME 

6. 15+00 

SECONDS 

POTENTIAL 

-1.15+04 

\'OLTS 

CYCLE 

23 

TIME 

6. 44+00 

SECONDS 

POTENTIAL 

-1. 15+04 

VOLTS 

CYCLE 

£4 

TIME 

6.73+00 

SECONDS 

POTENTIAL 

-1. 15+04 

VQLTS 

CYCLE 

85 

TIME 

7. 03+00 

SECONDS 

POTENTIAL 

-1. 15+04 

VOLTS 

CYCLE 

26 

TIME 

7.32+00 

SECONDS 

POTENTIAL 

-1.15+04 

'■'OLTS 


INCIDENT ELECTRON CURRENT -1.00-05 

SECONDARY ELECTRONS 8.19-06 

EPCfcSCATTERED ELECTRONS 1.81-06 
INCIDENT PPDTON CURRENT 0.00 

SECONDARY ELECTPONS 0.00 


NET CURRENT -4,06-10 AMPS.'M++2 


CYCLE 

27 

TIME 

7 . 61 + 00 

SECONDS 

POTENTIAL 

-1. 

CYCLE 

28 

TIME 

7. 90+00 

SECONDS 

POTENTIAL 

-1. 

CYCLE 

2 Q 

TIME 

8.20+00 

SECONDS 

POTENTIAL 

-1. 

CYCLE 

30 

TIME 

8.49+00 

■SECONDS 

POTENTIAL 

-1. 

CYCLE 

31 

TIME 

3.73+00 

SECONDS 

POTENTIAL 

-1. 


INCIDENT ELECTRON CUPPENT -1. 00-05 

SECONDARY ELECTRONS 9.19-06 

BACK SCATTERED ELECTRONS 1. 81-06 
INCIDENT PROTON CURRENT 0.00 

SECONDARY ELECTRONS 0.00 


15+04 VOLTS 
15+04 VOLTS 
15+04 VOLTS 
15+04 VOLTS 
15+04 VOLTS 


NET CUPPENT -1.11-11 AMPS/M++2 
PRINT A- CHARGING TABLE" Vi'ES 


T (SEC) 

0 . 00 
5.86-01 
1. 17+00 
- 1. 76+00 
2 . 34+00 
2.93+00 
3.51+00 
4. 1 0+00 
4.68+00 
5. 27+00 
5.86+00 
6.44+00 
7. 03+00 
7.61+00 
8 . 20+00 
8.78+00 
PLOT '' YES 


V (.VOLTS) 
0 . 00 

-2.34+03 
-4.58+03 
—6 . 66+03 
-8.51+03 

_o. .39+ 03 

-1.09+04 
-1 . 13+04 
-1 . 14+04 
-1. 14+04 
-1.15+04 
-I . 15+04 
-1 . 15+04 
-1.15+04 
-1.15+04 
-1.15+04 


I •.AMPS^'M++2) 

-7.26-06 
-7. 01-06 
—6 . 67—06 
-6. 16—06 
-5. 38-06 
-4. 14-06 
—2. 37-06 
-9.52-07 
-2.65-07 
-6.42-08 
-1.52-08 
-3.57-09 
-8.36-10 
-1.97-10 
—4. 62—1 1 
-1.11-11 


2. 00+04+ 

I 
I 
I 
I 
I 

1.50+04+ 

I 
I 
I 
I 
I 

1. 00+04+ 

I 
I 
I 

I ♦ 

I ♦ 

5. 00+03+ 

I + 

I ♦ 

I + 

I 

T ♦ 

~ J Q.OO 1.50+00 3. 00+00 4.50+00 6.00+00 7.50+00 9.00+00 

POTENTIAL ('-‘OLTS) VS TIME (SECONDS) 

EMISSION FORMULATION IS NORMAL 
SECONDARY EMISSION FORMULATION) 

MATERIAL 'STOP 

■EXIT) 

> 
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•3XQT 


WELCOME TO MATCHG - hi MATERIAL CHARGING PROGRAM 

EMISSION FORMULATION IS ANGLE 
SECONDARY EMISSION FOPMULATION> 

MATERIAL' KAPTDN 
PRINT > 

ENTEP FLUX TYPE 1=BEAM £=MRXI.JELLI AN 3 1 

BEAM VOLTAGE = 2. 0+01 KEV BERM CURRENT = 1.0-05 AMPS ✓ M++2' 

NEU PARAMETERS ? > 

GENERRTE R TRBLE 
CHARGE y >YES 
FULL PRINT r> 

CYCLE 1 TIME 0.00 SECONDS POTENTIAL 0.00 VOLTS 
INCIDENT ELECTRON CURRENT -1.00-05 

SECONDARY ELECTRONS 1.88-06 

BRCKSCRTTEPED ELECTRONS 5.70-07 
INCIDENT PPOTON CURRENT 0.00 

SECONDARY ELECTRONS 0.00 


NET CURRENT -7.55-06 RMP1/M++2 

CYCLE 31 TIME 2.46+03 SECONDS POTENTIAL -1.80+04 VOLTS 
INCIDENT ELECTPON CURRENT -1.00-05 

SECONDARY ELECTPONS 3.63-06 

BACK SCATTERED ELECTRONS 1.37-06 
INCIDENT PROTON CURRENT 0.00 

SECONDARY ELECTRONS 0.00 


NET CURRENT -3.71-12 AMPS/M++2 
PRINT A CHARGING TABLET > 

PLOT * s YES 


8. 00+04+ 
I 
I 
I 
I 
I 

1.50+04+ 

I 

I 

I 

I 

I 

1 . 00+04+ 


I * 

I 

I + 

I 

I + 

5. 00+03+ 

I + 

I 

I 

I +• 

I 

0.00 5.00+02 1,00+03 

POTENTIAL * VOLTS' 1 vs TIME 
EMISSION FORMULATION IS ANGLE 
SECONDARY EMISSION FORMULATION- 
MATERIAL. 

PRINT >YES 


♦ ♦ ♦ ♦■ ♦* ♦ ♦ ♦ +*• ♦ + ■*• 


1.50+ 03 

2. 00+03 

2.50+03 

3, 00+03 

(SECONDS; 
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PPEPROCESSIMG OF MATERIAL PROPERTIES 
MATERIAL It KfiPTOh 


PROPERTY 

DIELECTRIC CONSTANT 
THICKNESS 
CONDUCTIVITY 
ATOMIC NUMBER 
DELTA MAX >COEFF 
E-MBX >DEPTH++— 1 

RANGE 

EXPONENT > PENCE 

RANGE > EXPONENT 
EXPONENT 

YIELD FOP 1KEV PPQTONS 
MfiX DE'DX FOR PROTONS 
PHOTOCURPENT 
SURFACE PESISTIVITY 
CHANCE ANY MATERIAL PPOPEPTIES> 

ENTER FLUX TYPE 1=BEAM 2=MAXUELLIAN 
TE = 1 . O+OO KE« NE = 1 . 0+06 M-3 TI = 

NEW PPRAMETEPS T ^NO 
GENERATE A TABLE NO 
CHARGE t > YES 
FULL PRINT T>NO 


l 


4 

5 

6 
7 
S 

Q 

in 

H 

15 

13 

14 


INPUT VALUE 
3.50+00 CNONE) 

1.00- 04 METERS 

1.00- 14 MHO-' M 
5. 00+00 
5. 10+00 
1.50-01 

-1*00+00 ANG. 
0.0 0 (NONE) 
1.42+00 ANG. 
9.30+00 'NONE'* 
1.40+00 
7. 00+01 

2.00- 05 ft'M*+£ 
1.40+01 OHMS 


t,NONE> 

i.NaNE;i 

KEV 


<NONE> 
KEV 


>2 

1. 0+00 KEV NI 


CODE VALUE 
3.50+00 <*NQNE‘> 

1.00- 04 MESH 

1.00- 14 MHO^'M 
5 . 00+00 <NONE> 
3.05+01 -NONE) 
4.62-02 ANG- 01 
7.73+02 ANG. 

0. 00 ANG. 
1.51+00 <NONE> 

1. Cl 0+0 0 (HONE} 
1.40+00 CNCINE^ 
7. 00+01 KEV 

2.00- 05 A^M**2 
1. 24-10 V-S-cQ 


= 1 . 0+06 M-3 


CYCLE 


1 TIME 0. 00 


SECONDS POTENTIAL 0. 00 VOLTS 
T ilCI DENT ELECTPON CUPPENT -8.46-07 

SECONDARY ELECTPONS 1 . 24-06 

BACKSCATTEPEB ELECTPONS 2.32-07 
INCIDENT PROTON CURRENT 1.97-08 

SECONDARY ELECTPONS 

NET CURPENT 7. 15-07 AMPS^M^+2 

FLUX IS POSITIVE - WANT TO STOP? >NO llm „ 

CYCLE 3i TIME 2.60+00 SECONDS POTENTIAL 1.48+00 VOLTS 
INCIDENT ELECTRON CUPPENT -8.47-07 

SECONDAPY ELECTRONS 5.93-07 

BACKSCATTERED ELECTRONS 2.33-07 
INCIDENT PROTON CURPENT 1.97-08 

SECONDAPY ELECTPONS 3.40-08 


NET CURRENT 

PRINT A CHARGING TABLE 7 SYES 


3.20-08 AMPS^'M++£ 


V < VOLTS'* I »AMPS-'M++89 


0. 00 

ft. 00 

7. 15-07 

1.73-01 

3. 17-01 

5.23-07 

3.46-01 

5.55-01 

3.98-07 

5.20-01 

7.40-01 

3. 10-07 

6.93-01 

8.86-01 

2.47-07 

8 . 66—0 1 

1. 00+00 < 

! . 99-07 

1 . 04+00 

1 . 1 0+00 

1.62-07 

1.21+00 

1 . 18+0 0 

1.33-07 

1.39+00 

1.24+00 

1 . 1 0-07 

1.56+00 

1.30+00 

9.11-08 

1.73+00 

1.34+00 

7.60— OS 

1.91+00 

1.38+00 

6.37-08 

2. 08+00 

1 .41+00 

5.34-08 

2. 25+00 

1.43+00 

4.50-08 

2.43+00 

1. 46+00 

3.79-08 

2.60+00 
PLOT 7 >YE2 

1.48+00 

3.20-03 


**&%? 
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2 . 00 + 00 + 

I 
I 
I 
I 
I 

1. 50+00+- 
I 
I 

r 

i 
i 

1 . QO+00+ 

I 
I 
I 
I 
I 

5. 00— 01+- 
I 

I * 

I 

I + 

I 

rj ’"o.OO* 5. 00-01 1. 00+00 1.50+00 2.00+00 

POTENTIAL O'OLTS;- TIME '.SECONDS! 

EMISSION FORMULATION IS ANGLE 
SECONDARY EMISSION FORMULATIDN> 

MATEPIAL> 

PPINT 

ENTER FLUX TYPE 1=BEAM 2=MAXWELLIAN 
TE - 1.0+00 KEY ME = 1.0+06 M-3 TI = 1.0+00 LEV MI = 

MEM PARAMETERS T > 

GENERATE A TAELE t> 

CHARGE ? > YES 

FULL PPINT 

CYCLE 1 TIME 0.00 SECONDS POTENTIAL 0. UO 
INCIDENT ELECTPON CURPENT -8.46-07 

SECONDARY ELECTRONS 1.24-06 

BACKSCATTEPED ELECTRONS 2.32-07 
INCIDENT PPOTOH CIJPPENT 1.97-08 

SECONDARY ELECTPONS 7.11-08 




2.50+00 


1.0+06 M-"3 


VOLTS 


NET CURRENT 7.15-07 AMPS'M++2 
FLUX IS POSITIVE - MhNT TO STOP? >YES 
EMISSION FORMULATION IS ANGLE 
SECONDARY EMISSION FOPMULAT I ON> NORMAL 
EMISSION FORMULATION IS NORMAL 
MATERIAL! 

PPINT > 

ENTER FLUX TYPE 1=BEAM 8=M AXMELL I AN > 

TE = 1.0+00 KEY NE = 1 10+06 M-3 TI = 1.0+00 KEV NI = 1.0+06 M-3 

NEW PARAMETERS s 
GENERATE A TABLE 
CHARGE ? >YES 
FULL PPINT 7> 

CYCLE 1 TIME 0.00 SECONDS POTENTIAL 0.00 »OLTS 
INCIDENT ELECTRON CURRENT -3.46-07 

SECONDARY ELECTPONS 6.44-07 

BACKSCATTERED ELECTPONS £. 32-07 
INCIDENT PROTON CURRENT 1.97-08 

SECONDARY ELECTRONS 3.56-08 


3. 00+00 
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MET CUPPENT 3.58-08 HMPS.-M++2 


FLUX IS POSITIVE - MfitiT TO STORY .'YES 
EM I SSI DM FORMULATION I ? MORMPL 
-SECDNDAPY EMISSION FOPMULATIDNS 
MATERIAL.^ 

PRINT ^ 

ENTER FLUX TYPE 1=BEAM 2=MAXWELLIAN 
TE = 1.0+00 KEY ME = 1.0+06 M— 3 TI = 

NEW PPPBMETEPS ? > YES 

TE '.KEV> >5 
NE CM*+-3'> > 

TI <KEV> >5 
m |, M++— 3" 1 'i 

TE = 5.0+00 KE' J NE = 1.0+06 M-3 TI = 

GENEPBTE ft TABLE Y> 

CHftPGE y >YES 


l'. 0+00 KEV NI = 


5.0+00 KEV NI = 


FULL PRINT y> 

CYCLE 1 TIME 0.00 SECONDS POTENTIfiL 
INCIDENT ELECTRON CURRENT -1.89-06 

.SECONDARY ELECTRONS 6.40-0? 

BftCKSCftTTERED ELECTRONS 4.24-0? 
INCIDENT PROTON CURRENT 4.41-08 

SECONDARY ELECTRONS 1.58-0? 


0 . 00 


1. 0+06 M-3 


1. 0+06 M-3 


•'□LTS 


NET CURRENT -6.35-0? AMPS-'M*-^ 

CYCLE -31 TIME 7.43+03 SECONDS POTENTIfiL -3.68+03 VOLTS 
INCIDENT ELECTRON CURPENT -9.07-07 

SECONDARY ELECTRONS 3. 0?-0? 

BfiCKSCftTTEPED ELECTPGNS 2. 03-07 
INCIDENT PROTON CUPPENT 7.66-08 

SECONDARY ELECTPONS 2.98-n? 

NET CURRENT -2.22-08 AHPS.'M++2 
PRINT A rHARCINC TABLED )YES 


T <SEC> 

V *'VOLTS> 

I CftMPS/M++2> 

0. 00 

0.00 

-6.25-07 

4.95+02 

-.3. 06+02 

-4.68-0? 

9.90+02 

-1.42+03 

-3.59-07 

1 . 49-t-03 

-1.90+03 

-2.80-07 

1.98+03 

-2.27+03 

-2.21-07 

2.48+03 

—2 . 57 +03 

— 1 . 76—07 

£.97+03 

-2.31+03 

-1.41-0? 

3.47+03 

-3. 00+03 

-1. 14-07 

3.96+03 

-3.16+03 

-9.21-Oft 

4.46+03 

-3.28+03 

-7.43-08 

4.4S+03 

-3. 34+03 

-6. 08-08 

5.45+03 

-3.47+03 

-4 . 96-OS 

5.94+03 

-3.54+03 

-4. 05- OS 

6.44+03 

-3.59+03 

-3.31-08 

6.93+03 

-3.64+03 

-2.71-08 

7. 43^. ij-3 
PLOT S'* YES 

-3.68+03 

-2.22-08 



4 . 00 + 03 + 

I 

I 

I 

I 

I 

3 . 00 + 03 + 

I 

I 

I 

X 

i 

2 . 00 + 03 + 

I 

I 

I 

I 

I 

l . 00 + 03 + 

I 

I 

I 

I 

X 

m. on +- 

o. oo 


1.00+03 2.00+03 3.00+03 

PUTENTIAL <VQLT3'/ TIME 
EMISSION FOPMULRTIOH IS NORMAL 
SECONDARY EMISSION FDRM1JLAT ION> ANGLE 
EMISSION FORMULATION IS ANGLE 

material.^ 

PRINT 

ENTER FLUX TYPE 1=BERM 
TE = 5.0+00 KEV NE = 1. 

NEW PARAMETERS ? S 
GENERATE A TABLE '> 

CHARGE t > YES 
FULL PPINT TS 
CYCLE 1 TIME O.00 
INCIDENT ELECTPON CURRENT 

SECONDARY ELECTRONS 
3ACKSCATTEPED ELECTRONS 
INCIDENT PROTON CURRENT 


4.00+03 5.00+03 o. 00+03 7.00+03 t 

'.SECDNDS > 


2=MftXUELLIAN 
0+06 M— 3 TI = 


5.0+00 KEY NI = 1.0+06 M-3 


SECONDS POTENTIAL 
-1 . 39—06 
1.28-06 
4.24-07 
4.41-08 


0 . 00 


UOLTS 


SECONDARY ELECTRONS 


3. 16-07 


NET CURRENT 

FLUX IS POSITIVE - WANT tq 3T0R-? 

EMISSION FORMULATION IS ANGLE 
SECONDARY EMISSION FOPMULAT I ON> 
MATER I AL> STOP 


1.63-07 AMPS7M-++2 
-YES 


■'EXIT'J 


3FIN 


PUN ID: E48 ACCT: 11073-00 

TIME: TOTAL: 00:00:12.041 

CPU: 00:00:02.573 

CC-'EP: 00:00:08.059 

SUAS USED: 0.0419 

IMAGES READ: 141 PAGES: 

START: 09:14:14 JUL 13*1979 

♦TERMINAL INACTIVE* 


PROJECT: KATE— I 

CBSEC 

S: 00 0000323 

I-'O: 

00: 00: 01 .409 

WAIT: 

00:27:56.297 

£1 

FIN: 

09:42:29 JUL 13*197; 


!. 00+03 
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15. CDC 6600 CONVERSION 


Two versions of NASCAP were converted to run on the 
CDC 6600 machine. A version without booms, subdivision, emit- 
ters, or detectors was converted in April 1978. The up-to- 
date version was converted in October 1978. 

The task was completed in three stages. First, a 
machine independent set of NASCAP subroutines was created. 
Second, new versions of all machine dependent subroutines were 
created for the 6600. Third, the complete set of subroutines 
was loaded and executed on the 6600. 

15.1 MACHINE INDEPENDENT SUBROUTINES 

FORTRAN IV Extended Version on the CDC 6600 differs 
from FORTRAN V on UNIVAC machines in several ways. The most 
significant differences for the purpose of this conversion 
were : 

1. Non-executable Statements — order of these is 
ignored by UNIVAC, significant on CDC. Octal 
DATA statements also differ. 

2. PARAMETER Statements — used heavily in NASCAP, 
not allowed on CDC. 

3. Bit Manipulation — FLD in FORTRAN V, replaced by 
SHIFT and MASK on CDC. 

4. Data Transfer — different commands for efficient 
file communication. 

5. Graphics Routines — totally different graphics 
packages . 

6. Different machine word lengths. 

Our aim was to produce subroutine versions which were 
not affected by differences 1-6, for as many subroutines as 
possible. Differences 1 and 2 were easily handled in all 
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subroutines. We rearranged the order of statements and re- 
placed all parameter variables with their current values . 

Differences 3, 4, and 5 were handled by the introduc- 
tion of a new "layer" of subroutines. Bit manipulation (dif- 
ference 3) occurs in several dozen subroutines. Rather than 
make duplicate versions of these by replacing FLD with SHIFT 
and MASK, we created two new subroutines KBITS and SETBTS . 

We replaced each FLD bit manipulation with an equivalent call 
to KBITS or SETBTS. Now, any of the several dozen subroutines 
will run on either machines. We have duplicate versions of 
KBITS and SETBTS which use the UNIVAC or the CDC bit manipula- 
tion commands. Similarly, file input and output (difference 4) 
are handled by routines MOVDAT and CELLIO. All graphics rou- 
tines (difference 5) now refer to a set of seven basic sub- 
routines which are machine dependent. 

Difference 6 was more stubborn. In many places it 
seemed awkward to try to use the above technique. We had to 
settle for making some machine dependent subroutines by chang- 
ing " 10A6" formats to "6A10" formats. 

15.2 MACHINE DEPENDENT SUBROUTINES 

After granting independence to as many NASCAP sub- 
routines as possible, we were left with about forty machine 
dependent subroutines. Most of them fall into the following 
categories : 

1. "Extra layer" routines, as described in Section 
15.1. 

2. Routines with "6A10" changed to " 10A6" . 

3. Local system routines. 

4. A few routines with UNIVAC compiler commands. 

It was a straightforward process to rewrite these for 
the CDC 6600. 
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15.3 LOAD AND EXECUTE 


Loading of NASCAP on the 6600 required a new program 
segmentation to be devised. Execution naturally exposed a 
few oversights ,in the earlier phases of the process. How- 
ever,. both were accomplished without encountering particular 
problems. No unusual solutions were employed. 
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16. SCATHA CHARGING ANALYSIS 


Chapter 16 is a verbatim reproduction of a paper given 
at the 1978 USAF/NASA Spacecraft Charging Technology Conference. 
This paper is a clear summary of the SCATHA model analysis 
which was performed using NASCAP. 
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CHARGING ANALYSIS OF THE SCATHA SATELLITE* 


G. W. Schnuelle, D. E. Parks, I. Katz, 

M. J. Mandell, P. G. Steen, J. J. Cassidy 
Systems, Science and Software 

A. Rubin 

Air Force Geophysics Laboratory 
ABSTRACT 

We describe here a detailed model of the geometrical, 
material, and electrical properties of the SCATHA satellite for 
use with the NASA Charging Analyzer Program (NASCAP) . Charging 
calculations in an intense magnetospheric substorm environment 
demonstrate that: (1) long booms can significantly perturb the 

potentials near the spacecraft, and (2) discharging by sunlight 
or by active control can cause serious time-dependent differen- 
tial charging problems. 

1. INTRODUCTION 

We have developed a detailed model of the SCATHA satel- 

[ 1 . 2 ] 

lite for use with the NASA Charging Analyzer Program (NASCAP) . 

The model accounts for such geometrical complexities as booms , 
shadowing, and the presence of insulating materials over portions 
of the conducting ground of the space vehicle. The effects of 
photoemission and secondary emission caused by electron and ion 
impact, active control devices such as electron and ion beams, 
and surface and bulk conductivity are included in the model. To 
our knowledge, this model represents the most complete and re- 
alistic treatment of spacecraft charging attempted to date for 
any satellite. 

Section 2 below describes the SCATHA model employed in 
NASCAP. A detailed shadowing study was performed for a geo- 
metrically more accurate SCATHA model; this work is described 
in Section 3. We have performed charging calculations for one 

k _ _ 

This work supported by the National Aeronautics and Space Ad- 
ministration, Lewis Research Center, under Contract NAS3-21050. 
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environment using the present model, and the results of these 
calculations are described in Section 4. Preliminary conclu- 
sions of this study are summarized in Section 5. 

2. SCATHA MODEL DEVELOPMENT 

The NASCAP program allows the specification of the geo- 
metrical, material, and electrical properties of a spacecraft in 
considerable detail. We have attempted to incorporate the most 
current and complete information available for SCATHA into our 
model. However, the present model is meant primarily to illus- 
trate the intended level and scope of our study, rather than to 
provide the final word on a model specification. The NASCAP code 
allows model features to be easily altered to make our model a 
more faithful representation of the SCATHA satellite if the need 
arises . 

Perspective views of our gridded model are shown in Fig- 
ures 1 and 2. The main body of the satellite is represented as 
a right octagonal cylinder, with the aft cavity visible in Fig- 
ure 2. The OMNI antenna and the SC9 cluster of experiments are 
visible on the forward surface of the satellite. Our model re- 
produces the actual SCATHA geometrical features extremely well, 
as shown in Table 1. Note in particular that the treatment of 
booms in NASCAP allows the actual boom radii to be reproduced 
exactly in the model. The requirements in NASCAP that booms 
parallel coordinate axes and intercept mesh points in all grids 
effectively force any long booms to pass through the center of 
the innermost mesh. Therefore, our present model includes only 
the SC6, SC11, and the two SC2 booms, with the orientations 
fixed at right angles to one another. 

Figure 3 illustrates the computational space in which 
NASCAP solves Poisson's equation for this model. Monopole bound- 
ary conditions are imposed on the edges of the outermost grid, 
which is a rectangular prism of dimensions 1.6 x 1.6 x 3.2 m. 

The zone size decreases by a factor of 2 in each of the four 
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successive inner grids, so that the effective resolution is 
11.5 cm near the satellite body. (Local mesh refinement tech- 
niques in NASCAP allow a resolution of 2.5 cm for selected zones 
on the satellite.) 

Our model includes the specification of 15 distinct ex- 
posed surface materials, each of which is specified by the values 
of some 13 user-supplied parameters. The surface materials are 
described in Table 2 . We have attempted to find experimentally 
measured values for all parameters; where this has not been 
possible, suitable estimates based on the properties of similar 
materials have been used. Table 3 lists the values employed in 
the calculations reported here. The analytical expressions in 
which these parameters are used to evaluate net surface currents 
are described in detail in Reference 5. The formulation of 
electron backscattering in NASCAP has been somewhat modified 
recently, and the newer treatment is described in Appendix A. 

The exposed materials are illustrated in Figure 4 in which the 
locations of several of the SCATHA experiments are also shown. 
Experiments at the ends of SCATHA booms are modeled as a single 
boom segment whose radius is adjusted to match the exposed sur- 
face area of the actual experiment. 

The model includes six distinct underlying conductors: 
spacecraft ground, the reference band, and the four experiments 
SC2-1, SC2-2, SC6-1 and SC6-2. Each of these underlying con- 
ductors is capacitively coupled to spacecraft ground, and each 
can be separately biased with respect to ground. A seventh 
conductor could be introduced to underlay the solar cells at an 
appropriate bias. In this study the reference band was allowed 

to float and all other conductors were biased to the ground 

* 

potential. 

NASCAP has extensive capabilities to model particle emit- 
ters and detectors located on the spacecraft body, as described 
previously (Reference 2) . These features of NASCAP can be used 
in the analysis of the operation of, for example, the SCATHA 
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experiments SC4, SC5, SC6, SC7, and SC9 . Such studies should 
be particularly helpful in determining the influence of space- 
craft fields on particles emitted during active control, and 
in determining the source of particles seen at detector sites. 

3 . SHADOWING STUDY 

For the SCATHA shadowing study, we were required to gene- 
rate percent shadowing tables for various experiments. We were 
able to generate accurate tables using relatively small amounts 
of computer time: less than 5 minutes Univac 1100/81 time was 

required for a table of 7560 entries. 

Since the geometrical capabilities of the NASCAP shadow- 
ing routines are more general than the rest of the code, we were 
able to employ a SCATHA model for shadowing in which each experi- 
ment was treated geometrically in much finer detail than in the 
model described in Section 2. Figure 5 shows the level of de- 
tail in a perspective view of the ML12-7 experiments on the for- 
ward surface. Booms were placed at their actual locations on 
the satellite, and the 1 experiments at the boom ends were given 
a great deal of geometrical complexity. Figure 6 shows the SC2-1, 
SC1-4, and SC6-1 booms as they were resolved in the shadowing 
study . 

These detailed geometrical shapes were input to the usual 
NASCAP shadowing routines (HIDCEL) for table generation. The 
tables cover satellite rotation in 1° increments for the satel- 
lite plane deviations from the sun line of -5° to +5°. 

4. CHARGING CALCULATIONS 

The model was subjected to an extremely intense substorm 
described by a superposition of two Maxwellian plasmas with the 
following parameters : 


262 





10 0 eV 



The effects of ambient space charge were neglected in the 
solution of Poisson's equation here, since the mean satellite 
radius, r , is much smaller than the plasma Debye length, 

r„ ^ 100 cm 
s 

A^ * 700 J~ * 2200 cm 
D f n 

e 


r /A_ ^ 0.05 
s D 

There was no sunlight present in the first calculation de- 
scribed below. 

Potential contours during the initial overall charging 
phase ('V’lO ^ seconds) are shown in Figures 7 and 8. The ques- 
tion of whether booms have a significant effect on the sheath 
potentials is clearly answered by examining Figure 9 , which 
shows potential contours in a plane a half meter below the 
plane of the booms. Figure 10 shows similar contours in a 
calculation with the booms omitted; the distortion of con- 
tours by the booms is obvious. While the boom radii are 
small, ^2 cm, the effect on potentials is related to the boom 
capacitance, which varies only logarithmically with radius. 
This results in long range potential interactions from thin 
booms,, where the characteristic decay distance is closer to 
the boom length than to the boom radius. 

The rapid initial charging is followed by a much slower 
development of differential charging, as illustrated in Fig- 
ure 11. For this example the maximum differential developed 
after 22 seconds was 700 volts and the maximum field strength 


263 



in a dielectric layer was 24,000 volts/cm. Figure 12 shows 
contours in the plane of the booms after 22 seconds; note the 
differential charging developed at the boom ends due to varia- 
tions in the material properties between the experiments and 
the boom coatings.. 

The two Maxwellian description of the plasma leads to a 

low overall charging voltage of only -7 . 3 keV despite the 

presence of a plasma component with an electron temperature of 

40 keV. For the particular case we have studied here, low 

energy protons are being collected at an enormous rate and these, 

augmented by the secondary electrons they produce, balance the 

incident electron current. NASCAP uses a proton collection 

model in which the collection increases linearly with voltage, 

which is valid in the present case where r /A n is small, as dis- 

F41 5 u 

cussed by Laframboise. Table 4 shows the detailed current 

balance near equilibrium for the boom surface material in the 
presence of the double Maxwellian environment described above. 
Also shown in Table 4 is a similar breakdown for the same mate- 
rial subjected only to the high energy single Maxwellian com- 
ponent. The equilibrium potential is -32 keV in this case, in- 
dicating that the final potentials reached would have been much 
lower had we employed a single Maxwellian plasma model. For 
both plasma models, the final potentials reached will depend on 
the exact values employed for the proton and electron induced 
secondary yields. Great care should be exercised in the deter- 
mination of the values and associated error estimates for para- 
meters which affect the production of secondary electrons in 
these and similar calculations. 

Finally, the atomic number dependence of backscatter co- 
efficients tends to make high-Z materials charge less negatively 
than other elements. For SCATHA, this means that the magnitude 
of the boom potentials will be significantly lower than most 
other surfaces, since exposed platinum constitutes much of its 
surface area. 
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We have performed a similar calculation on this model 
in which the sunlight was turned on after 22 seconds of 
charging in eclipse. The photoemission results in strong 
differential charging (m3 keV) along the booms, as shown in 
Figure 13. In our model the boom surfaces are very weakly 
capacitively coupled to the grounded cable shields which ex- 
tend the length of the booms, while the experiments at the 
ends of the SC2 and SC6 booms are coupled closely to space- 
craft ground. This weak coupling has the effect of allowing 
the booms to react rapidly to environmental perturbations 
compared to the rest of the satellite, leading to temporary 
conditions of high differential charging. We have observed 
similar effects when discharging the satellite with an 
electron gun. 

The potentials near the satellite in sunlight are 
dominated by the monopole field of the spacecraft body. A 
photoemitting boom surface element can discharge only to the 
value of the local monopole potential, since further discharge 
is limited by immediate reflection of photoelectrons. This 
has the amazing consequence that the booms, strongly perturb- 
ing in eclipse, now seem to disappear in the potential con- 
tours near the satellite body. Note that significant dif- 
ferential charging in sunlight along the SC2 booms will cer- 
tainly persist at equilibrium due to large differences between 
the photoemission from surfaces on booms and on the SC2-1 and 
SC2-2 experiments. Our calculations neglect any effective 
surface conductivity parallel to the booms due to the 
presence of a photosheath. The surface conductivity features 
of NASCAP could easily be invoked to simulate this effect, 
which would reduce the magnitude of the differential charging 
observed here. 

The calculations reported here were performed on the 
UNIVAC 1100/81 computer at Systems, Science and Software. 

Each cycle of charging and solution of the potential equations 
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required approximately 15 minutes CPU time during differential 
charging, and 5 minutes CPU time when no differential charging 
occurred. Approximately 10 cycles of each type were required 
for the calculations reported here. We have developed a second 
SCATHA model for testing purposes in which the zone size is 
twice that of the model presented here and the booms are 
shortened; computer times are reduced by roughly 80 percent 
for this model, and all of the results described above can be 
observed in calculations using the smaller model. The half- 
scale model will be useful whenever fine resolution on the 
satellite surfaces is not required. 

5. CONCLUSIONS 

We have completed the development of a detailed model of 
the SCATHA satellite. Preliminary results from calculations in 
one magnetospheric environment indicate that: 

© The presence of a low energy component in a two 
Maxwellian description of the magnetospheric en- 
vironment reduces the maximum charging of a satel- 
lite relative to that found for a single Maxwellian. 

9 The booms have substantial impact on potentials 
near the spacecraft in eclipse. 

© The use of high atomic number coatings , such as 
platinum on the booms, may increase the severity 
of differential charging. 

© Discharging by sunlight or by active control may 
lead to transient increases in differential charg- 
ing along the booms due to the weak coupling of 
the booms to spacecraft ground. 

Our calculations demonstrate that the prediction of 
spacecraft potentials for SCATHA is an exceedingly complex 
problem, in which the full capabilities of the NASCAP treat- 
ment of geometrical features, material properties, and dynamic 
interaction with the environment are utilized. We plan to con- 
tinue this study of SCATHA using NASCAP with particular emphasis 
on boom perturbations and the effects of active control. 
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APPENDIX A.’ ELECTRON BACKSCATTER 

Electron backscatter is modeled in NASCAP as a function 

of electron energy and mean atomic number of backscattering 

fsi 

material. The formulation first used in NASCAP was valid 

only for low Z-materials. To remove this restriction we have 

f 6 1 

used a formula of Burke 1 to obtain the backscatter coeffi- 
cient for isotropically incident electrons as 

D 177 

n, = 0.475 Z * - 0..40 . (Al) 

The backscatter coefficient for normal incidence, ti , is then 
found by solving the equation 

ti 1 = 2 [1 - g o (l-£n ri o ) ] / (An n o ) 2 (A2) 

which comes from assuming the angular dependent backscatter 
, . [71 

coefficient to be 

n (0) = ti exp[-(An n Q ) (1 - cos0)] . (A3) 

[41 

The energy dependence is then taken to be 

n 0 (e) = Y(e) (n Q + 0.1 exp (-e/5)) (A4) 

( 0 e < 5 0 eV 

An (20 e)/An 20 50 eV < e < 1 keV 

1 e > 1 keV 

where e is in keV. 

The energy dependent r] Q from (A4) is then used in (A2) 
or (A3) to calculate the relevant backscatter coefficient. 
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Table 1. Comparison of Actual SCATHA Geometrical 
Features to Gridded MASCAP Model 

Zone Size = 4.53 in. (11.5 cm) 



SCATHA 

MODEL 

Radius 

33.6 inches 

32.0 inches 

Height 

68.7 

68.0 

Solar Array Height 

29 

27.2 

Bellyband Height 

11.3 

13.6 

SC9-1 Experiment 

9 . 2 x 6 x 8 

9.1 x 4.5 x 

SC6-1 Boom 

1.7 (radius) 

1.7 


118 (length) 

113.2 

Surface Area 

A 

2.16 x 10' sq. rn. 

2.11 x 10 4 

Solar Array Area 

o 

i — i 

X 

co 

• 

j-H 

1.15 x 10 4 

Forward Surface Area 

0.36 x 10 4 

0.34 x 10 4 
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Table 2. Exposed Surface Materials 


GOLD: 
SOLAR: 
WHITEN : 
SCREEN: 


YELOWC : 
GOLDPD : 

BLACKC : 
KAPTON : 
SI02 : 
TEFLON : 
INDOX : 
YGOLDC : 


ML12 : 


ALUM: 
BOOMAT : 


gold plate 

solar cells, coated fused silica 

non-conducting white paint (.STM K792) 

SC5 screen material, a conducting fictitious 
material which absorbs but does not emit charged 
particles 

conducting yellow paint 

88 percent gold plate with 12 percent conductive 
black paint (STM K7-48) in a polka dot pattern 

conductive black paint (STM K748) 

kapton 

Si 02 fabric 

teflon 

indium oxide 

conducting yellow paint (.50 percent) 
gold (50 percent) 

ML12-3 and ML12-4 surface, a fictitious, material 
whose properties are an average of the properties 
of the several materials on the ML12 surfaces 

aluminum plate 

platinum banded kapton 
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Table 3. 


Material Properties for Exposed Surfaces 


a 


b 





BLACKC 


- 

Property 

GOLD 

SOLAR 

WHITEN 

SCREEN 

YELOWC 

GOLDPD 

KAPTON 

1 

- 

4 . 00+00 

3.50+00 

— 

3.50+00 

— 

3.50+00 

2 

1.00-03 

1.79-04 

5.00-05 

1.00-03 

5.00-05 

1.00-03 

1.25-04 

3 

CO 

1.00-14 

5.90-14 

OO 

5.00-10 

OO 

1.00-14 

4 

7.90+01 

1.00+01 

5.00+00 

1.00+00 

5.00+00 

7.01+01 

5.00+00 

5 

8.80-01 

4.10+00 

2.10+00 

0.00 

2.10+00 

1.03+00 

2.10+00 

6 

8.00-01 

4.10-01 

1.50-01 

1.0 0+0 0 

1.50-01 

7.20-01 

1.50-01 

7 

8.30+01 

-1.00+00 

-1 .00+00 

1.00+01 

-1.00+00 

8 .30+01 

-1.00+00 

8 

1.63+00 

0.00 

0.00 

1.50+00 

0.00 

1 . 6 3 + 00 

0.00 

9 

3 . 46+01 

2.30+00 

1.05+00 

0.00 

1 . 0 5+00 

3.46 + 01 

1.42+00 

10 

7.00-01 

2.08+01 

9.80+00 

1.00+00 

9 .80+00 

7.00-01 

9.80+00 

11 

4.00-01 

1.36+00 

1.40+00 

0.00 

1.40+00 

4.00-01 

1.40+00 

12 

5.00+01 

4.00+01 

7.00+01 

1.00+00 

7.00+01 

5.00+01 

7.00+01 

13 

2.90-05 

2.00-05 

2.00-05 

0.00 

2.00-05 

2.90-05 

2.00-05 


SI02 

TEFLON 

INDOX 

YGOLDC 

ALUMIN 

BOOMAT C 

ML12 

1 

4 . 00 + 00 

2 . 00+00 


_ 

— 

2 . 00+00 

_ 

2 

2.75-04 

1.25-04 

1.00-03 

1.00-03 

1.00-03 

5.00-03 

1.00-03 

3 

2.75-12 

1.00-14 

oo 

oo 

00 

1.00-10 

CO 

4 

1.00+01 

1.00+01 

2.44+01 

4.20+01 

1 . 3 0 + 0 1 

6.34+01 

6 . 00+00 

5 

2.40+00 

3 . 00+00 

1.40+00 

1.49+00 

9.70-01 

1.86+00 

1.00+00 

6 

4.00-01 

3.00-01 

8.00-01 

4.80-01 

3.00-01 

5.90-01 

3.00-01 

7 

-1.00+00 

-1.00+00 

-1 . 00+00 

-1.00+00 

2.60+02 

8.30+01 

-1.00+00 

8 

0.00 

0.00 

0 . 00 

0.00 

1 . 30+0 0 

1.63+00 

0.00 

9 

1.02+00 

2.00+00 

7.18+00 

1.02+01 

2.40+02 

3.46+01 

2.00+00 

10 

2.00+01 

1.67+01 

5.55+01 

4.20+01 

1.73+00 

7.00-01 

1.20+01 

11 

1.40+00 

1.40+00 

1.36+00 

1.00+00 

1.36+00 

4.00-01 

1.40+00 

12 

7.00+01 

7.00+01 

4.00+01 

6.00+01 

4.00+01 

5.00+01 

7.00+01 

13 

2.00-05 

2.00-05 

3.20-05 

2.40-05 

4.00-05 

2.72-05 

2.10-05 




Table 3. (Continued) 


a The materials are described in Table 2. 

The thirteen properties are as follows (see Reference 4 and 
Appendix A for further details) » 

Property 1: Relative dielectric constant for in- 

sulators (dimensionless) . 

Property 2 : Thickness of dielectric film or vacuum 

gap (meters) . 

Property 3 : Electrical conductivity (mho/m) . The 

value 0 ° indicates a vacuum gap over a 
conducting surface. 

Property 4: Atomic number (dimensionless). 

Property 5 : Maximum secondary electron yield for 

electron impact at normal incidence 
(dimensionless) . 

Property 6 : Primary electron energy to produce 

maximum yield at normal incidence' (keV) . 

Properties 7-10: Range for incident electrons. Either: 

P P 

r o I 0 

Range = P ? E + P g E 

where the range is in angstroms and for 
the energy in keV, 

or 

P 7 = -1. to indicate use of an empirical 
range formula 

3 

Pg = density (g/cm ) 

P 10 ~ mean atoin:i - c weight (dimensionless) . 

Property 11: Secondary electron yield for normally 

incident 1 keV protons . 

Property 12 : Proton energy to produce maximum second- 

ary electron yield (keV) . 

Property 13 : Photoelectron yield for normally inci- 

dent sunlight (A/m2) . 

c The dielectric constant and thickness for the boom surfaces 
were chosen to reflect the effective capacitance to the under- 
lying cable shield. 
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Table 4. 


_5 2 

Components of Incident and Emitted Currents {10 A/m ) for Boom Surface 
Material Near Steady State. 


Double 

Maxwellian 


Potential 

-7000 Volts 

Incident Electrons 

-4.6 

Resulting Backscatter 

2.7 

Resulting Secondaries 

.7 

Incident Protons 

.6 

Resulting Secondaries 

.6 


Single 

Maxwellian 

-32,000 Volts 

-2.3 

1.4 ’ 

.4 

. 2 


.3 



SC6-1 



Figure 1. 


SCATHA model: side view. The 

SC1-4 boom are not included in 


50 m antenna and 
this model. 


the 
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Figure 3 . Computational space surrounding the SCATHA model , 
showing the nesting of the grids. The tic marks 
along the axes indicate the outer grid zone size; 
the zone size decreases by a factor of two in suc- 
cessive grids. 
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Figure 4b. SCATHA model with exposed surface materials il- 
lustrated. 
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SC1-1 
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HLl2 
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§§5 BCOMAT 
HH KAPTON 
S 102 
i SOLAR 


ML12-3.4.6 


Figure 4c. SCATHA model with exposed surface materials il- 
lustrated. 
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Figure 7 . Potential contours in a vertical plane through 
SCATHA center (only two of the four grids are 
plotted) . Note the contours extending into the 
aft cavity. Time ^10 — ^ seconds. Contours from 
-450 to -1250 volts in 50 volt steps. 
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Figure 8. Potential contours in a horizontal plane through 

SCATHA center. Time 'v-10"3 seconds. Contours from 
-300 to -1200 volts in 100 volt steps. The relative 
orientations of the booms is the same in later fig- 
ures. The dimples in the potential contours near 
the boom ends are artifacts associated with an im- 
perfect match of potential interpolation functions. 


C-H 
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Figure 9 


Potential contours in a horizontal plane 1 m below 
SCATHA center. Time ^10"^ seconds. Contours from 
-250 to -1150 volts in 50 volt steps. 
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Figure 10. Potential contours in a horizontal plane 1 m below 
SCATHA center for a model in which the booms have 
been removed. Time ^10“^ seconds. Contours from 
-300 to -1900 volts in 100 volt steps. 
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Figure 11. Spacecraft potential versus time for two points 
on SCATHA satellite. 
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Figure 12. Potential contours in a horizontal plane through 
SCATHA center, with differential charging along 
booms. Time v22 seconds. Contours from -2000 to 
-7000 volts in 500 volt steps. 
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Figure 13. Potential contours for sunlit case in a horizontal 
plane through SCATHA center. Time v38 seconds. 
Contours from -1000 to -7500 volts in 500 volt 
steps . 
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17. TWOD - TWO DIMENSIONAL SPACECRAFT CHARGING COMPUTER CODE 


17 . 1 INTRODUCTION 

Fo# the purpose of validating the ability of NASCAP to 
predict the presence of a sheath of photoelectrons and second- 
ary electrons. Systems, Science and Software has constructed 
a two-dimensional (R-0) spacecraft charging code capable of 
predicting equilibrium potentials and space charge densities 
about an infinitely long, dielectric-coated cylinder. The 
TWOD code features 

1. A flexible finite element formulation of the 
electrostatic potential problem. 

2. A fast ICCG (Incomplete Cholesky Conjugate 

(31 

Gradient) potential solver. 

3. Efficient particle tracking algorithms. 

4. Self-consistent formulations for space charge 
due to low energy emitted electrons. 

5. Linear (Debye length) screening for ambient plasma. 

6 . An approximate ef f ective-surf'ace-conductivity ■ 

treatment for charge transport in the photosheath; 
the conductivity value is determined consistent 
with the emitted current and external field. 

7. A first-order-implicit time-stepping treatment to 
promote stable convergence. 

The above features are either based on NASCAP ideas, extensions 
of techniques used in NASCAP, or formulations which might, at 
least in principle, be incorporated into NASCAP. In its 
final form, TWOD will treat material properties and environ- 
ment characteristics identically with NASCAP. 
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17.2 CODE STRUCTURE 

A block diagram of the preliminary version of TWOD is 
shown in Figure 17.1. This version does not include input or 
object definition sections, nor does it calculate currents 
based on material properties and environment specifications. 
Implementation of such features will be simple and straight- 
forward, in many cases involving direct transference of 
NASCAP routines. 

TWOD operates on a computational grid consisting of 
a central conductor surrounded by a ring of surface nodes and 
successive rings of space nodes. The space nodes are spaced 
uniformly in angle and are at arbitrary, specified radii, 
allowing fine resolution near the surface while including 
a large volume of space. 

TWOD operates in a time-stepping fashion, although 
procedures have been biased in favor of producing stable con- 
vergence toward equilibrium at the expense of accuracy in the 
time history. At each step, the Charge Section calculates 
incident and low-energy emitted currents to each sur face node. 
The l,ow-energy particles are tracked to find any barriers to 
their escape and their contribution to space charge. It is 
known that tracking of these particles to determine sur- 
face currents is an unstable procedure; therefore, following 

r 77 

a suggestion of Whipple an effective photosheath con- 
ductivity, calculated as described below, is used. Finally, 
the net current to each node is calculated, and an estimate 
is formed for the derivative of each node current with 
respect to its potential. These current derivatives serve 
to limit the voltage excursions, and are chosen to ensure 
stability. (See below.) 

Where space charge has an important effect, it is at 
best marginally convergent. Therefore, we use the ICCG 
potential solver to calculate the potential change due to 
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Initialization 
Form R-6 Grid 



Figure 17.1. Block diagram of TWOD code. 
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the change in space charge. The old space charge is then 
averaged with the newly calculated space charge so as to pro- 
duce a maximum space charge induced potential change of two 
volts . 

To complete the timestep we must solve 
C 

TWOD solves Eq. (17.1) "implicitly", i.e., with the integrand on 
the right hand side evaluated at the advanced time. Since the 
integrand can be very nonlinear, we approximate 

a(t) s; a (t^) (17.2a) 

J ± (t) = + jyv.(t) - V i (t 1 )J (17. 2b) 

where the coefficients of the diagonal matrix coefficients 
are estimated so as to ensure stability. (The space charge 
terms, Ap g , are calculated as described above. Debye shielding 
is included in the capacitance matrix C.) Equation (17.1) now 
becomes ■ - 


v(t 2 ) - 


V(ti> 


-f 


[J(t) + a (t) V(t) ] dt + 


A £s 


(17.1) 


[C - (J 1 + a) (t 2 - t x )] [V(t 2 ) - V(t x )] 

= [J(t 1 ) + aV(t 1 )] (t 2 - t-j^) + Ap s (17.3) 

Equation (17.3) is solved using the ICCG potential solver to 
obtain potentials for use in the next timestep. 
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17-3 EFFECTIVE PHOTOSHEATH CONDUCTIVITY 

In order to handle photosheath currents stably within 

the implicit potential solution it is necessary, following a 

[7 1 . 

suggestion of Whipple, to express these currents m the 

form 

J 0 (0) - a 0 (0 ) E 0 (8 ) . (17.4) 

Equation (17.4) is already a substantial, approximation, be- 
cause it is (a) linear, (b) local and (c) neglects "production 
gradient" transport. However, it is almost always the case 
that either a 0 (0) is large enough to maintain a nearly uni- 
form potential over a photoemitting area, or there is very 
little transport in the electron sheath. Therefore, Eq. (17.4) 
will always give the correct potential and space charge con- 
figuration, well within uncertainties in material properties 
and environment specifications. 

Figure 17.2 shows the geometry to be used in calculating 
photoconductivity. J 0 (Z) represents the current density of 
those electrons emitted to the right of the plane and landing 
to the left- We wish to calculate 

CO 0O 

/ J 0 (Z)dZ = f p(Z) v 0 (Z)dZ = P s v 0 (17.5) 

z o z o 

where p has dimensions of surface charge density. Neglecting 
s 

curvature effects, 

e . 

mm 

P s = f J(e z )t (£ z )de 2 (17.6) 

0 

where is the normal electron energy component, £ m ^ n is the 
minimum energy required for escape from the satellite, and 
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t(e„) is the flight time. The mean transverse velocity of 
L 

such a particle is 


V e z> 


t ( 


b f / 




(17.7) 


0 0 


so that 


e . 
mm 


/ V e)<iz - / 


J(e„) de„ - 
Z Z m 


/ d h / 


dt 2 E e (Z(t 2 )) 


(17.8) 


Formula (17-8) can be put into relatively simple form if we 
assume 


E = constant (>0) 

u 

E e = E e z~ 

We then obtain 


(17.9a) 

(17.9b) 


£ • ~ 
mm 


a e (E Z' Sitin' J/ Z o } 


/ 


e z J(e z )d£ z 


Z 0 


/ 


4 e, 


x / dx-L g ^ 


(17 . 10a) 


where 


2 - 1/2 

g(x,a) = (4a + a ) £n 


2 + ax t x V4a + a _ 
4 + 4ax(l - x) 


(17.10b) 
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£ „ is in electron volts, E„ is m volts/meter, and JCe,,) is 

“2 Z -2 . 
in amps/m -eV. Under most circumstances is the dominant 

factor in the conductivity. However, when E z becomes small, 

the integral also becomes small, so that the conductivity 

does not diverge. 
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17.4. EXAMPLE: ISOTROPIC FLUX 

As an example of the space charge feature of TWOD, we 

simulated a case in which the entire surface of a one meter 

radius cylinder was emitting photoelectrons at a rate of 
2 

1.25 nA/cm while receiving high energy incident electrons 

2 

at a rate of 0.125 nA/cm . Zero potential boundary conditions 
were set at a 12 meter radius. The equilibrium space charge 
and potential profiles are shown in Figure 17.3. The potential 
function clearly indicates formation of a space charge 
barrier about halfway between the surface and the outer 
boundary. The details of the potential profile are surely 
sensitive to the outer boundary condition, and would be modi- 
fied by Debye screening. (The space charge anomaly near the 
surface is due to poor statistics for the lowest energy parti- 
cles . ) 
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17.5 COMPARISON OF TWOD AND NASCAP 


It is important to recall that TWOD includes a self- 
consistent treatment of the space charge in the photosheath 
and an implicit effective surface conductivity arising from 
charge transport in the photosheath. NASCAP treats the ef- 
fects of photoemission implicitly by first estimating the 
maximum possible potential changes about a sunlit surface 
element and then adjusting the incident flux accordingly. 
Effective photosheath surface conductivity is not treated in 
NASCAP. The validity of these approximations is assessed 
by the comparisons given below. 

The TWOD program treats problems in R-0 geometry? we 

have performed calculations on a one meter radius cylinder 

with sunlight incident from one ‘side. The cylinder was 

-A 

covered with a 10 * m thick dielectric with s = 1, and zero 

boundary conditions were forced at a radius of 12 m. The 

incident current was from a plasma with n = n. = 3 cm 3 and 

e 1 2 

T e = T^ = 1 keV, and the photocurrent was 2 nA/cm for normal 
incidence. No secondary or backscattered electron currents 
were included. Figures 17.4 and 17.5 show the equilibrium 
potential distribution and space charge density. The shaded 
surface reached a potential of -2926 volts and the most sun- 
lit cell reached -1029 volts. An electrostatic barrier of 
approximately 1 volt formed about 15 cm above the most sunlit 
cell, leading to a maximum in the space charge density in the 
same region. 

The corresponding NASCAP calculation was performed on 
a right octagonal cylinder with an average radius of 1 m and 
a length of 3.24 m. The cylinder was covered with a dielectric 
1 cm thick with e = 1. The zone size was 27 cm, and monopole 
boundary conditions were forced at the edge of the second 
nested grid. Material parameters in NASCAP were adjusted to 
produce no secondary and negligible backscattered electron 


299 





Figure 17.5. Photosheath density from TWOD code. The maxi- 
mum space charge density is -5.2 x 10“H coul/m^. 
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current. The ambient plasma density was n g = n^ = 1 cm , 
and the plasma temperature was adjusted so that the cylinder 
would charge to -2930 volts in the dark; this led to T g = T^ = 
1.26 keV. The calculation was allowed to proceed until the 
potentials reached equilibrium, and the resulting potential 
distribution at the midpoint of the cylinder is shown in 
Figure 17.6. The back surface potential was -2950 volts, and 
the most sunlit node reached -970 volts. The NASCAP "SHEATH" 
option was then used to track low energy particles from the 
cylinder surface using the fixed equilibrium potential fields; 
plots of the resulting predicted photosheath space charge 
density are shown in Figure 17.7. Note the appearance of a 
maximum in the space charge density about one zone, or 27 cm, 
above the most sunlit cells. 

A comparison of Figures 17.4 and 17.5 with Figures 17.6 
and 17 . 7 indicate that the qualitative features of the NASCAP 
and TWOD predicted potential and space charge distributions 
agree well. A more detailed comparison is given in this sec- 
tion. 

The equilibrium surface potentials are plotted versus 
angle in Figure 17.8. The predicted differential charging 
was 1897 volts using TWOD and 1980 volts using NASCAP, an 
error of 6 percent. This discrepancy is due to the effec- 
tive photosheath surface conductivity, which was ignored 
in the NASCAP calculation. The oscillations in the dark 
side potentials from the NASCAP calculation are an arti- 
fact resulting from the procedures used to calculate average 
cell potentials for the charging algorithm and to share cell 
currents among nodes. We can eliminate these' oscillations 
by using revised averaging and sharing procedures with the 
current derivatives, (3J/3V) , as weighting factors. 

The radial dependence of the potentials are compared in 
Figure 17.9. Both calculations show a small potential bar- 
rier above the most sunlit cell; the location of the barrier 
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Figure 17.6. Potential contours (two grids) from NASCAP code. 

Contours from -500 to -2700 volts in 200 volt 
steps . 
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Figure 17.7. Photosheath density for inner grid from NASCAP 
"SHEATH" routine. Contours from 0.0 to -4.9 x 
10“H coul/m3 in steps of .62 x 10“H coul/m^ . 
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Figure 17.8. Equilibrium surface potential versus angle for 
a cylinder sunlit from one side. 
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Figure 17.9. 


R (meters) 

Potential versus radial distance for a cylinder 
sunlit from one side r in direction of incident 
sunlight. Cylinder radius = 1.00 meter. 


is nearer the surface in the TWOD results, at 15 cm, than in 

the NASCAP results, where it is nearer 30 cm. However, since 

the zone size in the NASCAP calculation was 27 cm, a barrier 

nearer than this distance cannot be resolved. The maximum in 

-11 3 

the space charge density was 5.2 x 10 coul/m using TWOD, 

-11 3 

and 5.1 x 10 coul/m using NASCAP, in good agreement. 

The above comparisons demonstrate that in conditions of 
strong differential charging, the approximations used in NASCAP 
to treat the effects of a photosheath introduce only small 
errors . 
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18 . THIN PLATES 


The thin plate algorithms were coded and tested in both 
potential and charging sections of the code. These routines 
are fully compatible with the LONGTIMESTEP option. Restric- 
tions on the thin plates include that surface cells on the 
underside not be subdivided. The present potential plotting 
algorithms do not include the potentials on the bottom of thin 
plates . 

18.1 -THIN PLATE EXAMPLES 

In order to test the algorithms several plate capacitor 
configurations were calculated. The results indicate that the 
accuracy of the representation is consistent with that of the 
other geometries treatable by NASCAP. For the special case of 
a parallel plate capacitor with a guard ring to minimize 
fringing fields, the results were accurate to within a few per- 
cent. Figures 18.1 through 18.4 show contours for several dif- 
ferent configurations. Figure 18.4 being the parallel plates 
with a guard ring. 
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potential contours along the x-z plaiic op y 
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Figure 18.1. Conducting parallel plates in a grounded tank. The left plot shows the 
two plates at the same potential; the right plot has them biased two 
hundred volts with respect to each other. 
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Figure 18.3. A differentially charged dielectric covered plate 
in sunlight. The extra contours on the photo- 
emitting surface are an artifact of the potential 
plot routines which do not presently have the 
capability to include potentials on the back side 
of plates. 
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Figure 18.4. Guard ring parallel plate capacitor. These are the same configuration as 
Figure 18.1, but the amount of charge on the central plate can be cal- 
culated separately. This allows accurate comparison to analytical theory 




18.2 DEFINITION OF THIN PLATES 


The definition of a thin plate is similar to that of a 
rectangular parallelepiped, with the exception that the ' SURFAC' 
designation is replaced by 'TOP' or 'BOTTOM'. Example: 


PLATE 



CORNER 


3 

DELTAS 


2 

TOP 

+Z 


BOTTOM 

-Z 


ENDOBJ 




1 4 

4 0 

TEFLON- 
KAPTON 


The distinguishing characteristic of a thin plate is 
that nodes interior to it are doubly defined. Thus, if an 
object contains more than one thin plate, a consistent speci- 
fication of "top" and "bottom" must be maintained. Also, be- 
cause of the way "bottom" points are treated in the potential 
solver, a volume cell touching the "bottom" of one thin plate 
may not touch the "top" of another. 
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